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Abstract 

This  research  has  been  aimed  at  developing  methods  to  predict  mechanical 
wear  of  sliding  bodies  at  high  velocities.  Specifically,  wear  of  test  sled  slippers  at 
the  Holloman  High  Speed  Test  Track  at  Holloman  AFB,  NM,  is  being  considered. 
Developing  a  numerical  model  to  represent  the  velocity  range  achieved  at  the  test 
track  is  infeasible,  so  numerical  modeling  techniques  must  be  adopted.  Previous 
research  has  made  use  of  finite  element  codes  to  simulate  the  high  velocity  sliding 
event.  However,  the  extreme  velocities  at  the  test  track  can  create  numerical  errors 
in  the  finite  element  codes.  To  avoid  the  numerical  errors,  an  Eulcrian-Lagrangian 
hydrocode  called  CTH  has  been  used  to  allow  for  a  velocity  range  of  200  to  1,500 
meters  per  second.  The  CTH  model  used  in  this  research  performs  plane  strain 
analysis  of  a  slipper  colliding  with  a  6  /irn  radius  semi-circular  surface  asperity. 

The  slipper-asperity  collision  event  creates  pressure  waves  in  the  slipper  which 
leads  to  failed  cells  and  worn  material.  Equations  have  been  derived  to  represent 
the  onset  of  plasticity  and  elastic  wave  speed  through  a  material  under  plane  strain 
conditions.  These  equations  were  validated  using  the  CTH  model.  Several  failure 
criteria  were  evaluated  as  possible  methods  to  estimate  damaged  material  from  the 
sliding  body.  The  Johnson  and  Cook  constitutive  model  was  selected  because  of  its 
ability  to  handle  high  strains,  strain  rates,  and  temperatures.  The  model  developed 
in  this  thesis  calculates  total  mechanical  wear  between  49.31%  and  80.87%  of  the 
experimental  wear  from  the  HHSTT  January  2008  test  mission. 
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The  Use  of  Various  Failure  Criteria 
As  Applied 
To  High  Speed  Wear 

I.  Introduction 

The  purpose  of  this  research  is  to  examine  the  interactions  of  sliding  bodies  at 
high  velocities  that  lead  to  wear.  In  order  to  accomplish  this,  the  wear  of  test  sled 
slippers  at  the  High  Speed  Test  Track,  (HHSTT),  at  Holloman  Air  Force  Base  (AFB) 
is  being  considered.  Principle  concepts  of  impact  and  wave  propagation  are  used  to 
evaluate  different  failure  criterion  within  the  material.  This  chapter  will  discuss  the 
goals  of  this  thesis,  as  well  as  provide  a  background  on  the  HHSTT.  Also,  previous 
research  in  the  formulation  of  wear  models  will  be  discussed. 

1.1  Objective  of  Research 

The  HHSTT  performs  a  variety  of  tests  at  high  velocities  using  a  rocket  sled 
system  that  rides  on  a  set  of  rails.  The  sled  is  attached  to  the  rails  using  slippers, 
which  are  described  in  greater  detail  in  Section  1.2.  The  rail  is  composed  of  AISI 
1080  steel,  whereas  the  slippers  are  made  of  VascoMax  300,  a  maraging  steel.  The 
HHSTT  engineers  would  like  to  be  able  to  estimate  the  amount  of  wear  of  each  slipper 
to  check  whether  it  will  reach  a  critical  thickness  before  the  end  of  the  test  run.  The 
goal  of  this  research  is  to  develop  a  numerical  model  to  predict  mechanical  wear  of  an 
HHSTT  slipper  as  it  slides  down  the  track.  These  numerical  models  take  advantage 
of  known  viscoplastic  characteristics  to  quantify  an  amount,  or  volume,  of  material 
that  has  reached  a  prescribed  failure  criterion. 

1.2  Holloman  High  Speed  Test  Track 

The  HHSTT  is  a  rocket  powered  sled  test  track  located  at  Holloman  AFB  in  New 
Mexico.  This  test  track  facility  is  used  for  a  variety  of  experiments  including  ground 
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level  aerodynamic  studies,  investigating  hypersonic  conditions,  munitions  testing,  and 
egress  systems.  In  April  2003,  a  land  speed  record  of  2,885  m/s  (6,453  miles  per  hour) 
was  set  at  the  facility. 

A  rocket  sled  train  used  for  these  high  speed  experiments  typically  consists  of 
several  pusher  stages  and  one  forebody  stage.  Rockets  are  attached  to  each  sled 
to  accelerate  the  train  to  the  desired  velocity.  The  forebody  sled  carries  the  test 
payload  and  required  instrumentation  along  with  the  final  rocket.  The  setup  shown 
in  Figure  1.1  is  the  configuration  used  for  a  mission  conducted  in  January  2008. 


Figure  1.1:  January  2008  Rocket  Test  Sled 


The  sleds  ride  on  parallel  AISI  1080  steel  rails  approximately  6,000  meters  long. 
Each  sled  is  attached  to  the  track  by  four  slippers  that  wrap  around  the  rail.  Figure  1.2 
shows  the  forebody  sled  attached  to  the  rail  for  the  record-setting  mission  in  April 
2003. 

Slipper  material  is  selected  based  on  the  maximum  velocity  of  each  slipper.  The 
first  two  pusher  sleds  in  the  January  2008  configuration  use  AISI  4130  steel  inserts 
placed  between  the  rail  and  slipper  housing,  as  shown  in  Figure  1.3.  The  steel  inserts 
are  discarded  after  each  test,  whereas  the  slipper  housings  are  reused.  The  third 
pusher  sled  and  forebody  sled  do  not  use  the  steel  inserts,  they  only  use  slippers 
fabricated  from  VascoMax  300.  These  slippers  are  discarded  after  each  test  due  to 
the  amount  of  wear.  Figure  1.4  shows  the  VascoMax  300  slippers  attached  to  the  rail. 
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Figure  1.2:  HHSTT  Rocket  Sled  System 


Figure  1.3:  VascoMax  300  Slipper  with  AISI  4130  Steel  Insert 
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Figure  1.4:  VascoMax  300  Slipper  without  Steel  Insert 
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Sled  designers  have  considered  several  variables  that  may  lead  to  poor  data 
collection,  or  failed  experiments.  One  concern  is  the  amount  of  wear  within  the 
slippers  as  they  reach  desired  velocities.  Figure  1.5  shows  the  standard  dimensions 
of  the  slippers  used  for  the  last  pusher  sled  and  forebody  sled.  The  designers  would 
like  to  be  able  to  estimate  the  amount  of  wear  of  each  slipper  to  check  whether  it  will 
reach  a  critical  thickness  before  the  end  of  the  test  run. 


Figure  1.5:  HHSTT  Slipper-Rail  Dimensions 

1.3  Summary  of  Previous  Research 

Research  into  the  mechanics  of  wear  has  been  ongoing  for  several  decades.  This 
research  has  led  to  several  definitions  of  wear.  The  American  Society  for  Testing  and 
Materials,  (ASTM),  defines  wear  as  “damage  to  a  solid  surface,  generally  involving 
progressive  loss  of  material,  due  to  relative  motion  between  that  surface  and  a  con¬ 
tacting  substance  or  substances”  [5].  Researchers  have  developed  various  methods 
and  experiments  to  define  the  rate  at  which  material  is  removed  as  it  slides  against 
another  surface. 

One  type  of  experiment  is  the  pin  on  disk  experiment.  This  experiment  uses  a 
rotating  disk,  or  ring,  and  pin  placed  on  the  surface  of  the  rotating  disk.  A  force  is 
applied  to  the  pin  and  the  material  removed  is  measured.  Figure  1.6  is  a  schematic 
of  a  pin  on  disk  experiment. 


5 


Figure  1.6:  Pin  On  Disk  Experiment  [13] 

In  1956,  Archard  and  Hirst  [4]  used  a  pin  on  disk  experiment  to  study  the  wear 
of  metals  under  unlubricated  conditions.  From  their  research,  Archard  and  Hirst 
concluded  that  wear  rate  was  initially  dependent  on  time,  until  interface  equilibrium 
was  reached  and  the  wear  rate  became  constant. 

In  1960,  Wolfson  [36]  studied  the  wear  of  materials  in  high  speed  track  appli¬ 
cations.  Sixty  tests  were  performed  at  varying  velocities,  bearing  pressures,  track 
conditions,  and  sliding  materials.  The  test  allowed  a  sled  to  accelerate  down  a  track. 
Once  the  desired  velocity  was  reached,  a  pin  was  dropped  into  contact  with  the  rail 
with  a  pneumatic  device.  This  pin  was  held  in  contact  with  the  rail  at  a  constant 
bearing  pressure  during  the  test.  The  pin  was  removed  from  the  rail  once  a  specified 
sliding  distance  was  reached.  The  amount  of  worn  material  was  determined  by  com¬ 
paring  final  dimensions  and  weights  to  initial  values.  Wear  rates  were  determined  by 
dividing  the  volume  of  worn  material  by  the  sliding  distance.  A  conversion  method 
was  applied  to  Wolfson’s  data  to  compare  results  to  the  analytical  model  developed 
in  this  thesis.  This  conversion  method  and  application  of  Wolfson’s  data  is  described 
in  greater  detail  in  Section  4.6. 

In  1970,  Farrell  and  Eyre  [19]  used  pin-on-disk  wear  experiments  to  characterize 
wear  between  two  steels.  Their  work  provided  distinction  between  mild  wear  and 
severe  wear,  the  transition  between  the  two  states,  and  its  dependence  on  both  sliding 
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speed  and  applied  load.  Mild  wear  “involves  the  relatively  slow  removal  of  the  tops  of 
the  highest  contacting  asperities  with  little  substrate  distortion,”  while  severe  wear 
shows  a  greater  scale  of  surface  damages  and  “the  wear  rate  increases  by  some  two 
orders  of  magnitude  from  that  of  mild  wear  and  the  maximum  size  of  the  wear  particles 
increases  suddenly  at  the  transition  load.”  Their  work  also  showed  that  the  coefficient 
of  friction,  /i,  is  dependent  on  both  the  sliding  velocity  and  applied  load,  as  shown  in 
Figure  1.7. 


Figure  1.7:  Variation  of  Coefficient  of  Friction  with  Load  [19] 


The  pin-on-disk  experiments  by  Archard  and  Hirst,  and  Farrell  and  Eyre  con¬ 
sidered  velocities  on  the  order  of  10  m/s.  In  1976,  Montgomery  [30]  published  the 
friction  and  wear  of  metals  in  high  muzzle  velocity  weapons.  A  pin-on-disk  experiment 
was  used  with  a  velocity  range  of  3  to  550  m/s.  In  order  to  keep  the  pin  from  running 
over  the  same  path  on  each  rotation  of  the  disk,  the  pin  was  moved  radially.  Strain 
gauges  were  used  to  measure  the  frictional  and  normal  forces  during  the  experiment. 
Similar  to  Farrell  and  Eyre,  Montgomery’s  experiments  showed  that  the  coefficient 
of  friction  was  dependent  on  both  the  sliding  velocity  and  applied  load,  or  pressure. 
The  coefficient  of  friction  was  plotted  as  a  function  of  the  product  of  pressure  and 
velocity,  Pv,  in  Figure  1.8. 

At  low  Pv  values,  the  coefficient  of  friction  was  higher.  As  the  Pv  value  in¬ 
creased,  the  coefficient  of  friction  decreased  exponentially  to  an  asymptotic  value. 
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Figure  1.8:  Variation  of  Coefficient  of  Friction  with  Pv  [30] 


According  to  Montgomery,  “the  mechanism  of  wear  at  high  sliding  speeds  is  almost 
certainly  surface  melting  followed  by  subsequent  removal  of  a  portion  of  the  melted 
surface  layer.”  This  surface  melting  creates  a  film  of  molten  material  along  the  sliding 
interface,  and  effectively  lowers  the  coefficient  of  friction  to  the  asymptotic  value  in 
Figure  1.8. 

In  1987,  Lim  and  Ashby  [26]  published  a  paper  describing  the  various  mecha¬ 
nisms  of  wear.  This  paper  formulated  wear-mechanism  maps  showing  the  relationship 
between  wear  mechanisms  and  test  conditions,  sliding  velocity  and  pressure.  These 
wear  mechanism  maps  were  generated  by  applying  two  converging  approaches.  The 
first  approach  was  to  plot  experimental  results,  and  identify  the  mechanisms  by  ob¬ 
servation.  The  second  approach  was  to  use  numerical  equations  describing  each  mech¬ 
anism.  The  two  methods  generate  a  map  showing  the  total  wear  rate  and  define  the 
contribution  of  each  wear  mechanism.  Figure  1.9  shows  the  wear-mechanism  map  for 
steel. 

Contours  of  constant  normalized  wear  rates  were  superimposed  on  fields  show¬ 
ing  the  regimes  of  dominance  of  different  wear  mechanisms.  There  were  discontinu¬ 
ities  in  the  contours  when  they  cross  the  held  boundaries  into  the  regimes  of  severe- 
oxidational  wear  and  melt  wear.  The  wear  rates  given  in  parentheses  were  the  values 
when  mild  wear  takes  place.  The  shaded  regions  indicated  a  transition  between  mild 
and  severe  wear  [26].  The  parameters  were  normalized  to  allow  specimens  of  various 
sizes  and  shapes.  Section  2.2  describes  the  normalization  in  greater  detail.  Very  little 
work  has  been  presented  in  the  past  that  stresses  the  relationship  between  wear  and 
wave  mechanics,  which  is  a  goal  of  this  research. 

1.4  AFIT  and  HHSTT  Wear  Research 

In  2007,  Cameron  [12, 13]  used  equations  developed  by  Archard,  and  Lim  and 
Ashby  to  characterize  the  wear  of  the  HHSTT  slipper  from  the  2003  test  run.  A  code 
was  written  to  analyze  dynamic  data  and  estimate  mechanical  and  melt  wear  depths. 
The  data  was  provided  by  Holloman  using  a  program  called  Dynamic  Analysis  and 
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Figure  1.9:  Wear  Mechanism  Map  for  Steel  [26] 
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Design  System,  (DADS).  The  DADS  data  is  discussed  in  greater  detail  in  Section  3.1. 
The  test  was  a  simulation  of  the  forebody  sled  accelerating  from  0  to  3,030  m/s  at 
a  constant  acceleration  for  2.5  seconds.  Cameron’s  analysis  calculated  a  mechanical 
wear  depth  of  0.27  cm  and  a  melt  wear  depth  of  0.08  cm.  The  total  wear  depth 
of  0.35  cm  is  less  than  the  nominal  thickness  of  the  HHSTT  slippers.  The  analysis 
was  deemed  an  acceptable  initial  approximation  of  high  velocity  slipper  wear  depth 
because  the  slippers  used  at  the  track  have  never  worn  through  the  entire  thickness 
due  to  a  test  run. 

In  2008,  Chrniel  [14]  used  a  finite  element  analysis,  (FEA),  approach  to  predict 
the  wear  of  HHSTT  slippers.  Two  methods  were  evaluated  in  the  research.  One 
method  used  equations  developed  by  Archard  on  a  macro-scale  in  incremental  steps, 
and  the  second  method  utilized  failure  criterion  based  on  material  property  on  a  micro¬ 
scale.  The  methods  were  evaluated  at  low  velocities  so  results  could  be  compared  to 
previously  published  experiments.  The  incremental  approach  produced  numerical 
errors  during  simulation  that  were  deemed  unacceptable.  The  failure  method  based 
on  material  properties  was  found  to  be  a  feasible  solution. 

In  2009,  Burton  [10, 11]  studied  the  surface  features  of  VascoMax  300  slipper 
and  AISI  1080  steel  rail  samples.  An  optical  profilometer  was  used  to  gain  accurate 
measurements  of  the  surface  roughness.  Figure  1.10  is  a  plot  of  recorded  surface  height 
data.  The  data  was  then  filtered  to  remove  surface  waviness  and  microstructural 
features.  Filtering  the  surface  data  was  beneficial  for  modeling  purposes  because  it 
removed  sharp  edges  and  sudden  changes  in  profile  which  can  lead  to  singularities 
when  used  in  FEA  models.  Figure  1.11  is  the  FEA  model  of  the  slipper  and  rail 
specimens  with  scanned  and  filtered  profile  geometry  used  by  Burton. 

The  FEA  model  was  used  to  study  the  effect  of  mesh  refinement  on  the  coefficient 
of  friction  at  the  interface  of  the  two  sliding  bodies.  Burton  stated:  “If  the  key  features 
of  the  surface  irregularity  are  not  represented  in  the  model,  the  model  coefficient  of 
friction  and  the  effective  coefficient  of  friction  for  the  macroscopic  forces  are  essentially 
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Figure  1.10:  VascoMax  300  Surface  Height  Data  [11] 


Figure  1.11:  Finite  Element  Model  Used  by  Burton  [11] 
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the  same.”  This  suggests  that  a  precise  measurement  of  the  surface  is  not  necessary 
to  model  the  coefficient  of  friction. 

In  2009,  Hale  [20]  used  a  micro-scale  FEA  approach  similar  to  Chmicl  to  model 
the  mechanical  wear  rates  of  a  hypothetical  HHSTT  test  run.  The  velocity  profile  of 
the  third  stage  from  the  January  2008  mission,  Figure  1.12,  was  used  for  this  research. 
The  wear  phenomenon  is  most  accurately  represented  as  a  3-dimensional  problem.  To 
simplify  the  model,  a  plane  strain  approach  was  used  to  simulate  a  VascoMax  300 
test  slipper  sliding  on  a  rail  made  of  AISI  1080  steel  and  colliding  with  a  semicircular 
surface  asperity  with  a  radius  of  6  /mi.  The  damage  criterion  used  was  based  upon 
the  viscoplastic  behavior  of  the  material  dehned  by  the  Johnson-Cook  [23]  model, 
discussed  in  Section  2.4.  The  total  damage  accumulated  by  each  finite  element  model 
(FEM)  run  was  divided  by  the  distance  slid  to  achieve  a  plane  strain  wear  rate. 


Third  Sled  Velocity  vs.  Time 


Figure  1.12:  HHSTT  Third  Stage  Velocity  Profile 


Hale’s  model  approximates  the  wear  of  HHSTT  slippers  by  a  collision  with  a 
single  surface  asperity  under  plane  strain  conditions.  In  a  real  test,  the  slipper  bounces 
and  slides  across  numerous  asperities.  In  order  to  account  for  the  multiple  asperities, 
a  scaling  factor  was  developed.  This  scaling  factor  was  determined  by  comparing 
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the  calculated  single  asperity  wear  rates  with  the  wear  rate  models  developed  by 
Archard  [2-4],  Applying  the  scaling  factor  allows  the  HHSTT  wear  problem  to  be 
simplified  to  a  simulation  with  a  single  asperity.  The  bouncing  of  the  slipper  was 
included  in  the  calculation  of  total  wear  by  multiplying  the  percentage  of  contact 
between  the  slipper  and  rail  during  a  test  run.  The  amount  of  slipper-rail  contact,  was 
determined  from  the  DADS  data.  Additionally,  a  coefficient  was  applied  to  represent 
a  semi-spherical  surface  asperity  in  the  plane  strain  simulation.  The  Archard  scaling 
factor,  contact  coefficient,  and  semi-spherical  coefficient  are  described  in  greater  detail 
in  Sections  3.3  and  3.4. 

In  2010,  Meador  [28]  used  a  hydrocode  to  investigate  the  wear  phenomenon. 
This  model  was  also  used  to  estimate  the  wear  of  a  hypothetical  HHSTT  test  run. 
However,  Meador  attempted  to  predict  slipper  wear  of  the  fourth  sled  reaching  a 
maximum  velocity  of  3,000  m/s.  For  this  research,  the  velocity  profile  was  identical 
to  the  third  stage  up  to  the  point  of  max  velocity  and  then  accelerates  to  3,000  m/s. 
Similar  to  Hale,  a  plane  strain  model  was  used  to  evaluate  failure  criterion  clue  to  the 
collision  with  a  surface  asperity. 

Meador  used  the  estimated  wear  rates  to  determine  the  total  wear  of  an  HHSTT 
slipper  for  an  entire  test  run.  The  total  wear  calculation  is  described  in  greater  detail 
in  Section  3.4.  The  results  of  this  calculation  were  compared  to  Hale’s  results  and 
experimental  data  from  the  2008  test  mission.  Figure  1.13  shows  the  total  wear 
volume  removed  for  a  sliding  distance  of  5,186  meters,  which  is  the  length  of  the 
January  2008  test  mission.  Meador’s  predicted  total  wear  was  greater  than  Hale’s, 
but  was  approximately  46%  of  the  total  measured  wear  from  the  2008  test  mission. 

In  2010,  Lodygowski  [27]  conducted  research  evaluating  the  temperature  of  two 
metals  sliding  relative  to  each  other.  The  FEA  model,  shown  in  Figure  1.14  forces 
a  plate  made  of  VascoMax  300  to  slide  between  two  AISI  1080  steel  surfaces.  The 
temperatures  of  the  material  were  calculated  over  a  velocity  range  from  1  m/s  to  200 
m/s. 
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Comparison  of  Mechanical  Wear  Calculations  with  Experimental  Wear 


Figure  1.13:  Total  Wear  Predicted  by  Meador  [31] 


Figure  1.14:  Finite  Element  Model  Used  by  Lodygowski  [27] 
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Lodygowski’s  research  provided  two  conclusions  relevant  to  this  thesis.  The  first 
conclusion  states  that  the  temperature  of  the  VascoMax  300  plate  does  not  change 
with  sliding  velocity.  This  is  due  to  the  fact  that  the  entire  VascoMax  300  plate  is 
not  in  contact  with  the  AIS1  1080  steel  throughout  the  whole  simulation.  However, 
the  interfacing  region  of  the  AISI  1080  steel  is  in  contact  with  the  VascoMax  300  for 
the  entire  simulation,  and  the  temperature  of  the  AISI  1080  steel  is  affected  by  the 
sliding  velocity.  The  model  developed  in  this  thesis  maintains  contact  between  the  two 
materials  during  the  entire  simulation.  As  such,  it  is  expected  that  the  temperature  of 
the  materials  will  be  affected  by  the  sliding  velocity.  The  second  conclusion  discusses 
the  relationship  between  material  temperature  at  the  interface  and  sliding  velocity. 
Figure  1.15  shows  the  average  temperature  of  the  AISI  1080  steel  along  the  interface 
for  a  given  sliding  velocity.  Lodygowski  states  that  the  relationship  is  not  linear,  but 
rather  a  logarithmic  increase  to  a  particular  value. 


Average  Temperatures  for  Specimen  A 


0  50  100  150  200  250 


Velocity  [m/s] 

Figure  1.15:  Variation  of  Temperature  with  Sliding  Velocity  [27] 


16 


1.5  Direction  For  Current  Research 


In  order  to  quantify  wear  of  HHSTT  slippers,  it  is  necessary  to  investigate 
the  mechanics  of  impact  that  lead  pressure  wave  propagation  through  the  material. 
There  has  been  much  consideration  of  impact  under  uniaxial  strain  conditions  and 
the  associated  pressure  waves  that  result  [29,37].  The  goal  here  is  to  extend  these 
understandings  to  plane  strain  scenarios  to  develop  a  model  that  represents  HHSTT 
environments,  but  to  stress  the  effect  of  wave  mechanics  in  an  associated  wear  envi¬ 
ronment.  This  is  discussed  further  in  Section  2.5.2. 

The  goal  of  this  research  is  to  create  a  model  that  will  accurately  predict 
mechancial  wear  of  VascoMax  300  slippers  colliding  with  a  of  6  /jrn  radius  semicircular 
surface  asperity  made  of  AISI  1080  steel.  Thus,  the  characteristics  of  wave  mechanics 
play  a  formidable  part  of  the  analysis.  A  hydrocode,  called  CTH,  discussed  in  Sec¬ 
tion  2.8,  is  used  to  simulate  this  scenario.  The  analysis  is  similar  to  the  micro-scale 
model  developed  by  Chrniel  [14],  and  used  by  both  Hale  [20]  and  Meador  [28].  Since 
a  numerical  model  is  used  to  evaluate  field  variables,  such  as  pressure,  stress,  strain 
rate,  etc.  there  are  several  failure  criteria  that  could  be  used.  Section  2.9  discusses 
the  various  failure  criteria  that  are  used  for  this  research. 
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II.  Theoretical  Background 


This  chapter  discusses  the  theoretical  background  required  to  develop  the  nu¬ 
merical  models  described  in  Chapter  III  and  interpret  the  results  presented  in  Chap¬ 
ter  IV.  A  description  of  the  various  wear  mechanisms  are  presented,  along  with  a 
discussion  of  the  coefficient  of  friction  between  sliding  metals.  This  chapter  will  also 
discuss  the  use  of  a  hydrocode,  including  the  considerations  of  conservation  equations, 
constitutive  equations,  and  an  equation  of  state.  Various  failure  criteria  used  to  quan¬ 
tify  material  damage  related  to  wear  is  presented  in  this  chapter.  Fundamentals  of 
wave  mechanics,  previously  defined  under  uniaxial  conditions  are  extended  for  appli¬ 
cation  in  a  plane  strain  scenario.  This  chapter  presents  derived  equations  defining 
the  onset  of  plasticity  and  elastic  wave  speed  through  a  material  under  plane  strain 
conditions. 

2.1  Wear  Rate 

The  model  developed  for  this  research  is  used  to  predict  mechanical  wear  rates 
defined  by  Equation  2.1  due  to  a  collision  with  a  surface  asperity,  where  W  is  the 
wear  rate,  Vw  is  the  volume  of  worn  material,  and  dsude  is  the  sliding  distance  into 
the  asperity.  Wear  rate  is  simply  dehned  as  the  volume  of  material  worn  per  distance 
slid.  Developing  a  model  to  predict  wear  rates,  allows  multiple  scenarios  to  be  run  and 
compared.  Specifically,  the  sliding  velocity  and  boundary  conditions  can  be  changed. 
It  was  found  by  Hale  [20]  that  the  wear  rate  from  a  mechanical  point  of  view  is  not 
history  oriented.  This  suggests  that  the  wear  rate  from  an  individual  simulation  is 
independent  of  wear  at  a  previous  simulation. 

W  =  ^~  (2.1) 

U'slide 

It  is  important  to  note  that  wear  is  a  system  response  influenced  by  both  mate¬ 
rial  properties  and  event  conditions.  These  event  conditions  consist  of  the  geometry 
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and  material  topography,  the  relative  motion  and  contact,  the  loading  scenario,  and 
any  environmental  conditions  including  lubrication  [6].  The  onset  of  wear  is  a  result 
of  collisions  between  surface  irregularities,  such  as  those  shown  in  Figure  2.1. 
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Figure  2.1:  Irregularities  in  Metal  Surface  Profile  [9] 


2.2  Wear  Mechanisms 

As  mentioned  earlier,  there  are  several  ways  to  define  wear.  In  order  to  avoid 
confusion,  for  the  purpose  of  this  research,  wear  is  defined  in  the  simplest  form  as  “the 
removal  of  material  volume  through  some  mechanical  process  between  two  surfaces” 
[31].  Furthermore,  there  are  several  mechanical  processes  that  can  lead  to  wear. 
Bayer  [8]  defines  three  ways  to  classify  wear.  In  no  significant  order,  the  first  is 
in  terms  of  the  appearance  of  wear.  The  surface  may  be  described  as  scratched, 
polished,  pitted,  etc.  The  second  classification  is  the  physical  mechanism  leading  to 
surface  damage.  Terms  related  to  this  classification  are  adhesion,  abrasion,  melting, 
and  oxidation.  The  third  classification  describes  the  situation  of  the  event  including, 
dry  sliding  wear,  lubricated  wear,  rolling  wear,  and  metal-to-metal  sliding  wear. 
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As  discussed  in  Section  1.3,  Lim  and  Ashby  [26]  developed  wear  mechanism  maps 
based  on  loading  scenario  and  material  properties.  The  wear  rates  were  normalized, 
W,  and  plotted  against  the  normalized  pressure,  F,  and  normalized  velocity,  v.  The 
wear  rate,  pressure,  and  velocity  were  normalized  using  Equations  (2.2,  2.3,  and  2.4) 
respectively. 

SLIDING  VELOCITY  v  (m/s) 


NORMALISED  VELOCITY  v 

Figure  2.2:  Wear  Mechanism  Map  [26] 
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In  these  equations,  W  represents  the  wear  rate,  An  represents  the  normal  contact 
area,  F  represents  the  applied  load,  H  represents  the  material  hardness,  v  represents 
the  sliding  velocity,  ro  represents  the  radius  of  the  pin  used  for  the  experiment,  and 
a  represents  the  thermal  diffusivity  of  the  material.  The  normalization  equations  are 
used  to  relate  to  studies  using  different  size  and  shape  specimens.  Research  into  the 
various  wear  mechanisms  has  yielded  two  scenarios  of  interest  for  this  thesis;  abrasive 
wear  and  adhesive  wear,  both  of  which  fall  into  the  classification  of  mechanical  wear. 
All  research  presented  herein  is  only  considering  the  phenomenon  of  mechanical  wear. 
Other  wear  mechanisms,  such  as  melt  wear  and  oxidation  are  not  considered  in  this 
thesis. 

2.2.1  Abrasive  Wear.  Abrasive  wear  occurs  when  asperities  along  the  inter¬ 
face  of  the  sliding  bodies  collide.  The  tangential  force  is  large  enough  to  cause  plastic 
deformation  and  eventually  remove  the  asperity.  Figure  2.3A  represents  an  abrasive 
wear  scenario.  Material  from  the  asperity  is  being  removed  by  the  triangular  shaped 
abrasive  particle. 

2.2.2  Adhesive  Wear.  Adhesive  wear  occurs  when  two  surfaces  contact  at 
an  asperity  and  bond  together.  As  the  sliding  motion  continues,  and  if  the  bond 
is  strong  enough,  asperities  from  the  softer  material  will  shear  off  and  adhere  to 
the  harder  material.  The  adhered  fragments  later  break  free  forming  worn  material. 
Figure  2.3B  depicts  the  adhesive  wear  event. 

2.3  Coefficient  of  Friction 

As  two  bodies  slide  relative  to  each  other,  they  are  inhibited  by  friction.  Friction 
is  a  phenomenon  resulting  from  tangential  motion  between  the  two  sliding  bodies, 
and  conventionally  is  thought  of  as  the  force  required  to  initiate  or  to  sustain  the 
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Figure  2.3:  Abrasive  Wear  and  Adhesive  Wear  [7] 

tangential  motion  [28].  It  is  important  to  consider  the  frictional  forces  between  the 
rail  and  slipper,  as  wear  is  dominated  by  the  interface  of  the  two  materials.  The 
wear  mechanisms;  abrasive  wear,  and  adhesive  wear,  discussed  in  Section  2.2,  are 
relative  to  the  friction  in  the  sliding  bodies.  Furthermore,  the  local  temperatures  of 
the  materials  are  affected  by  frictional  heating  as  the  sliding  motion  occurs. 

Establishing  a  coefficient  of  friction  is  a  common  way  to  represent  the  frictional 
forces  between  two  surfaces.  The  coefficient  of  friction,  /i,  relates  the  frictional  force 
to  the  normal  force  applied  between  the  two  surfaces.  Equation  2.5  is  used  to  solve  for 
the  coefficient  of  friction,  where  Fj  is  the  frictional  force  and  F  is  the  normal  contact 
force.  The  equation  assumes  that  the  coefficient  is  independent  of  the  contact  area 
and  proportional  to  the  normal  load. 


h  = 


h 

F 


(2.5) 
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The  assumptions  listed  for  Equation  2.5  are  generally  accepted  for  mild  slid¬ 
ing,  or  low  velocity  motion.  However,  the  contact  area  begins  to  make  a  significant 
contribution  to  wear  as  the  sliding  velocity  and  loading  increases.  Research  by  Mont¬ 
gomery  [30],  presented  in  Section  1.3,  discusses  the  relationship  between  coefficient 
of  friction  and  Pv.  Experimental  results  show  at  low  Pv  values,  the  coefficient  of 
friction  was  higher  for  steel-on-steel  sliding.  As  the  Pv  value  increased,  the  coefficient 
of  friction  decreased  exponentially  to  an  asymptotic  value.  Hale  [20]  applied  a  curve 
fit  to  the  tabulated  data  from  Montgomery  for  steel  sliding  on  steel.  This  curve  is 
used  to  represent  the  coefficient  of  friction  for  the  VascoMax  300  slipper  sliding  along 
the  AISI  1080  steel  rail  as  a  function  of  the  Pv  term.  Figure  2.4  shows  the  data  and 
curve  fit.  Equation  2.6  is  the  exponential  curve  fit.  It  is  important  to  note  that  the 
curve  fit  was  generated  using  Pv  data  with  units  of  MPa  •  nnn/s.  Any  use  of  the 
curve  fit  requires  the  same  units. 


Figure  2.4:  Montgomery  Data  with  Curve  Fit  [20] 
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2-4  Johnson- Cook  Viscoplasticity  Model 

In  1983,  Johnson  and  Cook  [23]  studied  metals  subjected  to  large  strains,  high 
strain  rates,  and  high  temperatures.  Test  data  for  the  model  was  obtained  using 
torsion  tests  and  dynamic  Hopkinson  bar  tensile  tests  over  a  range  of  temperatures. 
The  elevated  temperatures  were  obtained  by  surrounding  the  specimen  with  an  oven 
for  several  minutes  prior  to  testing.  Adiabatic  heating  resulting  from  high  strains 
complicated  the  results  because  the  elevated  temperatures  showed  a  reduction  in  the 
material  strength.  Adiabatic  heating  occurs  when  the  pressure  of  a  material  increases 
due  to  the  motion  of  surrounding  particles.  In  this  case,  high  strains  caused  the 
temperature  of  the  material  to  increase  without  adding  heat.  Johnson  and  Cook 
developed  Equation  2.7,  a  constitutive  model  to  solve  for  the  flow  stress,  a. 


a  =  [A  +  Be”]  [1  +  C  ln(e*)]  [1  -  T*m]  (2.7) 

This  equation  is  a  product  of  three  terms.  The  Erst  term  is  the  static  yield 
strength  and  a  modification  for  strain.  The  second  term  introduces  strain  rate  de¬ 
pendency  and  the  final  term  includes  temperature  effects  [20].  A,  B ,  C,  m,  and  n  are 
material  constants,  ep  is  the  equivalent  plastic  strain,  e*  is  the  dimensionless  plastic 
strain  rate  for  s^1. 

Equation  2.7  is  used  to  represent  yielding  from  the  effective  stress,  von  Mises 
stress.  The  stress  indicates  the  yield  surface  when  it  reaches  the  Johnson- Cook  equa¬ 
tion,  Equation  2.7.  At  that  point  a  corresponding  plastic  strain  rate  is  determined. 
Since  the  process  is  time  integrated,  the  incremental  step  time  of  a  simulation,  At,  is 
incorporated  in  the  calculation  for  the  incremental  strain.  This  incremental  step  time 
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is  discussed  in  Section  3.2.  The  Bodner-Partom  relationship  is  then  used  to  update 
stress  [35].  This  method  allows  subsequent  stress-strain  relations  to  be  developed. 
The  homologous  temperature,  T*,  is  defined  by  Equation  2.8. 


T*  =  T  ^  (2.8) 

Where  T  is  the  material  temperature,  T0  is  the  ambient  temperature,  and  Tmeit 
is  the  material  melting  temperature.  The  homologous  temperature  must  be  defined 
in  order  to  create  strain  rate  dependent  stress-strain  curves.  Due  to  the  high  strain 
rate  deformation  applied  for  this  research,  the  deformation  work  is  considered  adi¬ 
abatic.  This  implies  that  the  deformation  work  is  transformed  into  heat  with  the 
rise  in  temperature  of  the  material.  This  temperature  rise  is  observed  in  stress-strain 
curves  as  thermal  softening,  a  behavior  which  constitutive  equations  must  account 
for.  Meyers  [29]  defines  the  adiabatic  temperature  rise  in  a  material  subjected  to  high 
plastic  strain  rate  due  to  plastic  strain  energy  as  Equation  2.9. 

8  f 

AT  =  ——  /  ode  (2.9) 

pCP  Jo 

where  (3  is  the  inelastic  heat  fraction,  p  is  density,  Cp  is  the  specific  heat  ratio,  and 
ep  is  the  final  plastic  strain.  The  inelastic  heat  fraction  is  set  as  0.9  in  the  analysis 
based  on  results  from  ductile  materials  [29].  Equation  2.10  is  given  by  replacing  the 
stress  term  in  Equation  2.9  with  the  Johnson-Cook  constitutive  equation,  Equation 
2.7,  and  assuming  the  strain  rate  is  constant. 


dT*  _  /3(1  +  Cln(e;))  y6/ 
1  ~T*m  ~  pCp(Tmelt  -  Tref)  J0 


(A  +  B(ep)n)d£p 


(2.10) 


where  T0*  and  Tf  are  the  initial  and  final  homologous  temperatures.  Even  though 
m  is  a  material  constant,  m  =  0.8  for  VascoMax  300  [15,16],  Meyers  proposes  an 
approximation  of  m  ~  1  to  the  left  hand  integral  in  Equation  2.10.  Applying  Meyers’ 
approximation,  the  homologous  temperature  reduces  to  Equation  2.11. 
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(2.11) 


rj~ i* 


1  —  exp 


/i(l  +  Cln(ffl»))  /  B(eT^\ ' 
pCp(Tme[t  -  Tref)  V  n  +  lj 


The  true  stress-strain  curves  can  now  be  generated  by  substituting  Equation 
2.11  for  the  homologous  temperature  in  Equation  2.7.  Figure  2.5  shows  the  true 
stress-strain  curves  for  VascoMax  300  with  increasing  strain  rates.  The  Johnson  and 
Cook  constitutive  model  uses  variables  commonly  found  in  computational  software 
which  makes  it  easy  to  use  for  simulation. 


Figure  2.5:  True  Stress-Strain  Curves  for  VascoMax  300  with 
Johnson-Cook  Constitutive  Equation  [20] 


2.5  Wave  Propagation 

The  research  presented  in  this  thesis  is  an  investigation  in  the  wear  of  VascoMax 
300  due  to  the  collision  with  a  surface  asperity  made  of  AISI  1080  steel.  The  goal 
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is  to  characterize  mechanical  wear  at  high  velocities.  Analysis  of  a  high  velocity 
collision  requires  some  background  information  about  wave  mechanics,  specifically 
the  formulation  of  a  shock  wave  within  a  solid  medium.  The  next  few  sub-sections 
present  important  information  pertaining  to  the  mechanics  of  wave  propagation. 

2.5.1  Uniaxial  Strain.  Common  analysis  of  wave  propagation  through  solid 
media  has  been  studied  while  considering  uniaxial  strain  [29].  The  yield  point  for 
uniaxial  strain  is  referred  to  as  the  Hugoniot  Elastic  Limit,  written  as  cthel •  This 
is  the  maximum  elastic  stress  for  one-dimensional  elastic  wave  propagation  in  plate 
geometries  [37].  The  Hugoniot  Elastic  Limit  represents  the  onset  of  plasticity  in 
a  material  under  uniaxial  conditions.  Due  to  the  strain  limitations,  the  onset  of 
plasticity  is  greater  than  the  yield  stress,  Y0)  dehned  for  uniaxial  stress  conditions. 
The  Hugoniot  Elastic  Limit  is  given  by  Equation  2.12,  where  v  is  the  poisson  ratio. 

°HEL  =  Fo  (brt)  (2'12) 

It  is  important  to  determine  the  Hugoniot  Elastic  Limit,  because  if  the  stress 
in  the  material  exceeds  this  limit,  two  waves  are  created.  First,  an  elastic  wave  will 
move  through  the  material  at  the  clastic  speed,  Ce,  defined  by  Equation  2.13,  where 
E  is  the  clastic  modulus,  and  p0  is  the  initial  density  of  the  material.  Following  the 
clastic  wave,  a  plastic  wave  will  move  through  the  material  at  the  plastic  wave  speed, 
Cp,  speed  defined  by  Equation  2.14.  It  is  important  to  note  that  cp  is  a  function  of 
the  slope  of  the  stress-strain  curve  at  a  given  point.  This  means  that  multiple  waves 
can  exist  in  the  material,  each  with  a  speed  defined  by  Equation  2.14. 


ce 


E{l-u) 


Po(  1  -  2z/)(l  +  v) 


(2.13) 


Cp 


(2.14) 
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Figure  2.6  shows  the  transition  of  a  pressure  wave  to  a  shock  wave  with  in¬ 
creasing  time.  Point  A  represents  a  low  pressure  region  in  the  wave  moving  at  a  low 
velocity.  Points  B  and  C  have  higher  pressures  and  therefore  move  at  a  greater  veloc¬ 
ity.  The  four  steps  in  the  figure  show  the  pressure  wave  approaching,  and  ultimately 
reaching,  a  vertical  line  which  represents  the  onset  of  a  shock  within  the  material.  Be¬ 
fore  the  shock  is  formed,  the  material  properties  across  the  pressure  wave  are  smooth 
and  easily  defined.  As  the  shock  wave  propagates,  the  material  movement  becomes 
discontinuous  in  front  of,  and  behind  the  shock.  An  equation  of  state,  EOS,  is  used 
to  estimate  the  pressure  and  internal  energy  of  the  material.  Equations  of  state  are 
discussed  in  greater  detail  in  Section  2.8.2. 


Figure  2.6:  Increasing  Pressure  Wave  with  Time  [37] 

2.5.2  2D  Plane  Strain.  The  analysis  in  this  thesis  considers  plane  strain 
conditions.  The  approach  used  to  characterize  pressure  waves  under  uniaxial  condi¬ 
tions  is  modified  to  include  strain  in  two  dimensions.  In  making  these  modifications, 
an  assumption  was  made  that  the  strain  component  can  be  represented  by  the  sum¬ 
mation  of  the  two  principal  strains.  The  equivalent  Hugoniot  Elastic  Limit  for  plane 
strain  conditions,  (Thel,pSi  is  given  by  Equation  2.15.  This  equation  represents  the 
onset  of  plasticity  in  a  material  under  plane  strain  conditions.  Appendix  A  includes 
the  full  derivation  of  the  equations  presented  in  this  chapter. 


VHEL.PS  -  Y0  +  1  (2.15) 

This  scenario  does  restrict  strain  in  the  third  principal  axis  (z-direction) ,  while 
allowing  strain  in  the  other  dimensions.  Therefore,  the  value  of  <Jhel,ps  should  be 
less  than  <Jhel  but  greater  than  the  yield  stress,  Ya.  Table  2.1  shows  the  uniaxial 
yield  stress,  Hngoniot  clastic  limit,  and  equivalent  plane  strain  Hugoniot  elastic  limit 
for  VascoMax  300  with  a  Poisson’s  ratio  ( v  =  0.283).  The  values  for  yield  stress  and 
poisson  ratio  for  VascoMax  300  are  taken  from  Cinnamon  [15]. 

Table  2.1:  VascoMax  300  Hugoniot  Limits 


Y0  (GPa) 

&hel  (GPa) 

&hel,ps  (GPa) 

2.1 

3.4692 

2.8664 

The  speed  of  an  clastic  wave  under  the  plane  strain  condition,  ce,ps  is  given  by 
Equation  2.16.  The  speed  of  the  plastic  wave  under  plane  strain  conditions  is  still 
defined  by  the  slope  of  the  stress-strain  curve  at  a  point.  Multiple  plastic  waves  are 
still  produced  in  plane  strain,  each  defined  by  Equation  2.14. 


2.6  FEA  and  Hydrocodes 

Recent  research  at  AF1T  has  relied  on  the  use  of  two  codes  to  model  the  wear  of 
HHSTT  slippers.  One  is  an  FEA  code  called  ABAQUS  and  the  other  is  a  hydrocode 
called  CTH.  Both  codes  can  be  used  to  create  2-dimensional  or  3-dimensional  geome¬ 
try,  to  represent  a  wide  variety  of  simulations,  and  both  codes  make  use  of  a  mesh  to 
solve  the  numerical  analysis.  There  is  one  fundamental  difference  between  the  FEA 
approach  and  the  hydrocode  approach.  The  difference  stems  from  the  meshes  used 
and  frames  of  reference  established  in  each  code. 
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The  FEA  method  uses  a  Lagrangian  mesh  which  attaches  the  mesh  to  the  mate¬ 
rial.  This  means  the  mesh  will  deform  with  the  material  during  the  analysis,  and  the 
frame  of  reference  moves  with  each  successive  iteration.  This  method  is  an  attractive 
approach  for  many  simulations  because  the  equations  are  simple  to  solve,  ffowever, 
if  the  scenario  involves  large  deformations  leading  to  excessive  material  displacement, 
the  FEA  method  begins  to  fall  apart.  With  large  deformations,  numerical  singulari¬ 
ties  in  the  finite  element  equations  can  result  due  to  the  mesh  cell  geometry.  In  some 
cases,  the  element  can  invert  under  large  distortions,  resulting  in  negative  volume  and 
negative  mass  [37].  This  gives  rise  to  numerical  errors. 

Hydrocodes  use  an  Eulerian  mesh  which  fixes  the  mesh  in  free  space  and  allows 
the  material  to  flow  through  it.  The  frame  of  reference  does  not  move  during  the 
analysis.  Figure  2.7  is  a  simple  depiction  of  the  slipper-rail  scenario  with  a  Lagrangian 
mesh  (left)  and  an  Eulerian  mesh  (right).  The  Eulerian  mesh  is  commonly  referred 
to  as  a  finite  area  mesh. 


Lagrangian  Mesh 


Eulerian  Mesh 


Figure  2.7:  Graphical  Comparisson  of  Lagrangian  and  Eule¬ 

rian  Meshes  [28] 
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2.1  ABAQUS 


Previous  research  at  AFIT  used  the  Lagrangian  FEA  approach  to  model  a 
VascoMax  300  slipper  colliding  with  a  AISI  1080  steel  surface  asperity.  Initial  work 
by  Chmiel  [14]  used  ABAQUS  to  evaluate  Lagrangian  codes  as  a  way  to  model  the 
wear  of  HHSTT  slippers.  The  slipper  was  modeled  using  4-node  plane  strain  elements 
and  the  rail  used  3-node  plane  strain  elements.  Chmiel’s  model  is  shown  in  Figure  2.8. 
The  analysis  showed  material  damage  as  a  result  of  the  collision  with  the  asperity. 
This  proved  that  the  FEA  method  could  be  used  to  predict  mechanical  wear  of  the 
HHSTT  slippers.  However,  it  has  been  deemed  impractical  to  run  the  simulation 
the  entire  length  of  the  steel  rails  at  Holloman  AFB.  Chmiel  suggested  simulations 
of  single  asperity  collisions  over  a  range  of  velocities,  and  calculating  the  total  wear 
from  these  runs. 


Figure  2.8:  Finite  Element  Model  Used  by  Chmiel  [14] 


Hale  [20]  continued  the  research  of  HHSTT  slipper  wear  using  ABAQUS  to 
implement  Chmiel’s  suggestion.  The  asperity  collision  event  is  best  described  as  a 
3-dimensional  event,  so  a  3-dimensional  model  was  considered.  However,  complica- 
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tions  arose  when  evaluating  a  3-D  model.  A  plane  strain  model  was  used  to  calculate 
the  wear.  A  3-dimensional  conversion  method  and  scaling  factor  based  on  equations 
developed  by  Archard  were  used  to  extend  the  plane  strain  results  to  a  3-dimensional 
value.  The  3-dimensional  conversion  method  and  Archard  scaling  factor  are  described 
in  greater  detail  in  Section  3.4.  Hale’s  model,  shown  in  Figure  2.9,  used  a  combina¬ 
tion  of  3-node  linear  plane  strain  triangular  elements  and  4-node  bilinear  reduced 
integration  elements  to  model  the  slipper  and  rail. 


Figure  2.9:  Finite  Element  Model  Used  by  Hale  [20] 


2.8  CTH 

CTH  is  a  Lagrangian-Eulerian  hydrocode  developed  by  Sandia  National  Lab¬ 
oratories.  Hydrocodes  can  be  very  useful  when  evaluating  scenarios  involving  high 
velocity  impact  resulting  in  wave  propagation  and  the  possiblity  of  a  shock.  The  pro¬ 
cess  in  solving  a  hydrocode  solution  is  less  straight-forward  than  the  FEA  method. 
However,  there  are  detailed  descriptions  of  CTH  features  published  by  Sandia  Na¬ 
tional  Laboratories  [18]  and  Palazotto  and  Meador  [31].  Additionally,  Zukas  [37]  and 
Meyers  [29]  each  wrote  books  providing  descriptions  of  hydrocodes.  The  following  sec- 
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tions  use  information  taken  from  these  sources  to  describe  in  detail  the  key  features 
of  the  code,  and  how  it  is  used  to  model  the  HHSTT  wear  problem. 

2.8.1  Lagrangian  Step  and  Eulerian  Remap.  The  conservation  equations  of 
mass,  momentum,  and  energy  are  satisfied  using  a  two  step  process  in  CTH  [18,28,31]. 
First,  the  Lagrangian  step,  is  used  to  evaluate  the  equations  across  the  time  step  and 
the  mesh  deforms  with  the  material.  This  is  followed  by  an  Eulerian  remapping 
step  which  redefines  the  mesh  to  the  original  Eulerian  coordinates.  Equations  (2.17, 
2.18,  and  2.19)  are  the  Lagrangian  conservation  equations  for  mass,  momentum,  and 
energy,  respectively,  where  p  is  the  material  density,  is  velocity,  P  is  pressure,  a  is 
stress,  E  is  energy,  and  Q  represents  additive  heat  as  a  function  of  velocity  and  wave 
speed,  cs  [18].  These  equations  are  a  by-product  of  the  Eulerian  expression  in  which 
the  substantial  derivative  is  applied. 

f  (2-17) 

P ^  =  -  VP  -  V  •  [a  +  Q(rf,  cs)]  (2.18) 

P(l^  =  -PV  •  -  [a  +  Q(lf,  cs)}  •  W  (2.19) 

Since  the  mesh  deforms  initially,  the  conservation  of  mass  is  trivially  satisfied, 
because  no  mass  flow  occurs  across  the  cell  boundaries.  The  momentum  and  energy 
integrals  are  solved  using  their  explicit  finite  volume  representations  [28,31].  Thermal 
energy  of  the  material  must  be  considered.  The  conservation  of  energy  equation 
includes  both  mechanical  and  thermal  energies,  and  thus  another  equation  relationship 
that  couples  the  energies  together  is  needed.  This  is  where  an  equation  of  state  is 
used.  The  conservation  equations  and  equation  of  state  are  solved  in  conjunction 
with  a  constitutive  model.  CTH  decomposes  the  total  stress  tensor  into  a  spherical 
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part,  solved  for  using  an  equation  of  state,  and  a  deviatoric  part,  solved  for  using  a 
constitutive  model. 

After  the  conservation  equations  have  been  satisfied  and  the  constitutive  equa¬ 
tion  has  been  applied,  the  Eulerian  remap  step  is  used  to  return  the  distorted  mesh 
to  the  original  Eulerian  mesh.  An  interface  tracking  algorithm  internal  to  CTH  is 
utilized  to  track  locations  of  material  interfaces  within  mixed  cells  containing  mul¬ 
tiple  materials.  The  change  in  mass  is  calculated  by  the  geometry  of  the  deformed 
material  compared  to  the  previous  step.  The  mass  and  internal  energy  are  mapped  to 
the  fixed  mesh.  The  results  from  the  interface  tracking  algorithm  are  used  to  map  the 
momentum  and  kinetic  energies  to  the  material  in  the  Eulerian  mesh.  The  equation 
of  state  is  used  to  update  the  pressure,  temperature,  and  density  of  the  cells. 

2.8.2  Equation  of  State.  The  equation  of  state  is  used  to  relate  the  inter¬ 
nal  energy,  and  pressure,  of  a  material  to  the  density  and  temperature.  There  are 
several  equations  of  state  that  can  be  used  within  CTH.  For  the  purpose  of  model¬ 
ing  the  collision  of  two  solid  bodies,  CTH  provides  two  separate  EOS  models,  the 
semi-empirical  Mie-Griineisen  EOS,  and  the  tabular  Sesame  EOS.  Vanderhyde’s  [34] 
research  provides  insight  to  the  two  EOS  models  internal  to  CTH.  Much  of  the  infor¬ 
mation  presented  in  the  following  sections  use  information  taken  from  this  source. 

2.8.2. 1  Mie-Griineisen  EOS.  The  Mie-Griineisen  equation  of  state  is 
typically  used  for  high  velocities  ranging  from  500  m/s  to  2,000  m/s  [37].  No  phase 
change  is  allowed  with  the  Mie-Gruniesen  equation  of  state,  which  makes  it  useful  for 
this  research.  Equation  2.20  is  the  Mie-Griineisen  equation  used  by  CTH,  where  Pref 
and  Eref  are  reference  pressure  and  energy,  usually  taken  from  the  Hugoniot  relations 
or  by  assuming  a  zero-Kelvin  isotherm,  V  is  volume,  and  T  is  the  Griincisen  constant 
defined  by  Equation  2.21,  where  K  is  the  bulk  modulus,  v  is  the  specific  volume,  a  is 
the  thermal  coefficient  of  expansion,  and  Cv  is  the  specific  heat. 
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(2.20) 


P  Pref  —  y(E  Eref ) 


3  Kva 
Cv 


(2.21) 


2. 8. 2. 2  Sesame  EOS.  CTH  also  provides  the  Sesame  equation  of  state. 
The  Sesame  EOS  is  a  set  of  tabular  data  collected  through  experimentation  at  Sandia 
National  Laboratories  primarily  using  flyer  plate  impact  experiments  [15].  When  the 
Sesame  EOS  is  used,  CTH  interpolates  between  the  tabulated  data,  or  extrapolates 
outside  of  the  provided  data  set  to  estimate  the  internal  energy  and  pressure  within 
the  material. 


2.8.3  Boundary  Conditions.  Boundary  conditions  are  determined  using 
finite  volume  approximations  based  on  the  surrounding  cells.  However,  the  cells  along 
the  boundary  have  at  least  one  side  with  no  neighboring  cell.  Boundary  conditions 
need  to  be  established  in  order  to  solve  the  hnite  volume  problem  when  using  a 
hydrocode.  These  conditions  are  based  upon  the  concepts  of  sound  waves,  which  is  a 
primary  relationship  in  this  research,  and  thus  the  boundary  conditions  must  be  able 
to  control  mass,  momentum,  and  energy  fluxes  across  the  boundary.  There  are  four 
possible  boundary  conditions  in  CTH:  a  symmetrical  boundary  condition  (Type  0), 
a  sound  speed  based  absorbing  boundary  condition  (Type  1),  an  outflow  boundary 
condition  (Type  2),  and  an  extrapolation  boundary  condition  (Type  3).  The  boundary 
conditions  create  additional  cells  just  outside  the  internal  mesh  defined  in  the  problem 
setup. 

The  Type  0  boundary  condition  sets  parameters  of  the  adjacent  cells  equal  to 
the  cells  along  the  boundary  of  the  internal  mesh.  The  velocity  between  the  two 
adjacent  cells  is  set  to  zero  and  kinetic  energy  is  converted  to  internal  energy.  Also, 
mass  flux  is  restricted  across  the  boundary.  The  Type  1  boundary  condition  allows 
mass  to  enter  the  internal  mesh,  and  is  used  to  approximate  semi-infinite  bodies.  The 
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Type  2  boundary  condition  sets  empty  cells  on  the  boundary  of  the  internal  mesh 
with  user  specified  pressure.  Mass  can  exit  the  internal  mesh,  but  no  mass  is  allowed 
to  enter.  The  Type  3  boundary  condition  places  cells  on  the  boundary  of  the  internal 
mesh  and  linearly  extrapolates  a  boundary  pressure.  No  restrictions  on  mass  flux  are 
present  in  this  boundary  condition.  Previous  work  by  Meador  [28]  used  the  type  1 
boundary  condition.  The  analysis  presented  in  this  thesis  used  a  combination  of  type 
1  and  2. 

2.8.4  Data  Collection.  There  are  two  methods  of  recording  data  during 
a  CTH  simulation.  One  method  uses  locations  attached  to  the  stationary  mesh  to 
record  data  as  the  material  deforms  through  it.  The  other  method  uses  tracer  points 
that  travel  with  the  material  during  the  simulation.  For  the  purpose  of  this  research, 
the  second  method  is  used.  Utilizing  this  method  requires  the  use  of  a  tracer  input  set 
defined  within  the  CTH  input  deck.  The  tracer  input  set  defines  the  initial  locations  of 
each  data  point.  As  the  simulation  occurs,  the  tracer  points  move  with  the  material. 
This  method  keeps  track  of  failed  material,  representing  wear,  through  the  entire 
simulation.  If  the  data  were  collected  at  the  stationary  mesh  locations,  one  cell 
could  be  considered  damaged,  or  failed,  at  a  previous  time  step,  when  new  material 
has  entered  the  cell.  As  a  result,  this  new  material  will  not  be  qualified  as  damaged, 
because  that  cell  has  previously  been  defined  as  damaged.  This  results  in  unreasonably 
low  wear  predictions. 

2. 9  Failure  Criteria 

Properly  assessing  wear  requires  established  failure  criteria  to  quantify  material 
damage.  Meador  [28]  outlines  four  failure  criteria  to  determine  wear:  average  strain 
rate,  point-wise  strain  rate,  Johnson-Cook  plasticity,  and  plastic  strain.  The  Johnson- 
Cook  constitutive  equation  is  used  to  define  the  plastic  strain  failure  criterion,  and 
was  selected  to  evaluate  wear  in  the  model  developed  for  this  research.  A  second 
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failure  criterion  was  established  based  upon  a  critical  stress  value.  The  Johnson-Cook 
fracture  model  [24,33]  was  also  evaluated  as  a  failure  criterion. 

2.9.1  Plastic  Strain  at  Max  Stress  Failure  Criteria.  This  method  evaluates 
the  plastic  strain  at  maximum  stress  for  a  given  strain  rate  from  the  true  stress-strain 
curves,  as  given  in  Figure  2.5.  These  curves  were  developed  using  the  Johnson-Cook 
constitutive  equation,  Equation  2.7,  with  Meyer’s  approximation  for  the  homologous 
temperature  [29],  presented  in  Section  2.4.  Each  curve  on  the  plot  represents  the  true 
stress-strain  relationship  for  VascoMax  300  with  a  given  constant  strain  rate.  The 
critical  stress,  crcrit,  is  defined  as  the  maximum  stress  of  each  curve.  Similarly,  the 
critical  plastic  strain,  £prit ,  is  defined  as  the  strain  at  maximum  stress  for  each  curve. 
The  critical  strain  can  be  determined  as  a  function  of  the  strain  rate  by  plotting 
epcrit  against  the  strain  rate  and  applying  a  curve  fit  through  the  data.  The  curve  fit 
provides  a  closed  form  solution  for  the  critical  strain,  Equation  2.22. 


£crit(x,  y,  t )  =  APSe{x ,  y ,  t)Bps  +  CPS  (2.22) 


The  constants  AP$,  BPS,  CPs  are  given  in  Table  2.2.  Where  the  ‘PS’  subscript 
is  used  to  identify  plastic  strain  constants.  Figure  2.10  shows  the  plastic  strain  curve 
fit  for  VascoMax  300. 


Table  2.2:  Coefficients  of  Plastic  Strain 


Coefficient 

Value 

Units 

APs 

2.247  x  10"2 

MPa 

BPs 

-5.516  x  10"2 

unit  less 

Cps 

6.044  x  10”3 

MPa 

CTH  calculates  and  records  the  plastic  strain  at  each  tracer  point  during  the 
simulation.  A  MATLAB  post-processing  code,  written  to  compare  the  recorded  plastic 
strain  against  the  critical  plastic  strain  from  Equation  2.22,  is  discussed  in  greater 
detail  in  Appendix  C. 
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Critical  Plastic  Strain, 


Johnson-Cook  Plastic  Strain  at  Max  Flow  Stress  vs.  Plastic  Strain  Rate 


Figure  2.10:  Critical  Plastic  Strain  Curve  Fit  [28] 
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2.9.2  Von  Mises  Stress  Failure  Criteria.  CTH  can  be  used  to  calculate  the 
von  Mises  stress  of  the  material  during  the  simulation.  Therefore,  a  failure  criterion 
can  be  established  based  on  a  critical  stress  value.  This  method  would  allow  CTH  to 
do  all  the  calculations  in  determining  material  failure.  Hale  [20]  provides  maximum 
stress  values  for  VascoMax  300  based  on  a  dominant  strain  rate.  Table  2.3  shows 
the  dominant  strain  rates  and  associated  maximum  stress  from  the  Johnson- Cook 
constitutive  equation  for  a  range  of  sliding  velocities  from  the  January  2008  test 
mission.  Although  the  stress  levels  change  from  2,900  MPa  to  3,130  MPa,  the  critical 
stress  value  for  the  von  Mises  stress  failure  criterion  was  chosen  to  be  3,000  MPa. 
This  means  that  if  the  von  Mises  stress  exceeds  3,000  MPa  during  the  simulation,  it 
has  failed. 


Table  2.3:  VascoMax  300  Maximum  Stress  Based  on  Dominant  Strain 
[20] 


Velocity  Range  (m/s) 

Dominant  Strain  Rate 

Maximum  Stress  (MPa) 

10  -  200 

1  x  105 

2,900 

300  -  622 

1  x  106 

3,000 

750  -  1,530 

1  x  107 

3,130 

2.9.3  Johnson-Cook  Fracture  Model.  CTH  allows  the  use  of  the  Johnson- 
Cook  fracture  model  for  evaluating  failed  material.  “It  uses  a  failure  criterion  based 
on  equivalent  plastic  strain,  taking  into  account  the  pressure,  temperature,  and  strain 
rate  along  the  loading  path  for  each  material  particle.  The  model  uses  one  scalar 
damage  variable”  [33].  Equation  2.23  defines  the  plastic  strain  at  failure. 


/  -P3p  \ 

£pf(p,Y,T,e)  =  [D1  +  D{2  Y  ’Wl  FDA\n{max{\, e))][l  FDbT*]  (2.23) 

Where  p  is  pressure,  Y  is  the  material  yield  stress,  T  is  temperature,  e  is  the 
plastic  strain  rate,  D\ ,  D2,  D3,  D 4,  and  D5,  are  material  constants  derived  from  ex¬ 
perimentation,  and  T*  is  the  homologous  temperature  previously  defined  by  Equation 
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2.8.  The  Johnson-Cook  scalar  damage  variable,  D ,  is  defined  by  Equation  2.24.  Ini¬ 
tially,  undamaged  material  has  a  D  value  equal  to  0.  As  the  simulation  occurs,  the 
material  accumulates  damage,  and  the  scalar  variable,  D,  increases.  When  D  equals 
1,  the  material  is  damaged. 


D 


dep 

£pf(p,  Y,  T,  e) 


(2.24) 


2.10  Summary  of  Theoretical  Background 

The  information  presented  in  this  chapter  has  been  crucial  in  understanding  the 
wear  phenomenon.  The  relations  presented  will  be  used  to  develop  a  numerical  model 
to  predict  mechanical  wear  rates  of  HHSTT  slippers.  Previous  work  by  Hale  [20]  and 
Meador  [28]  has  made  use  of  a  plane  strain  scenario  to  model  the  slipper-rail  sliding 
event.  These  models  allow  a  VascoMax  300  slipper  to  collide  with  a  6  /irn  surface 
asperity  made  of  AISI  1080  steel.  Damage  was  recorded  per  sliding  distance  to  give 
wear  rates.  Wave  action  in  a  plane  strain  scenario  was  also  evaluated. 

The  fundamental  approaches  to  characterizing  the  onset  of  plasticity  in  the  uni¬ 
axial  strain  case,  (?hel,  were  presented.  Steps  were  taken  to  extend  the  characteristics 
to  the  case  of  plane  strain.  This  resulted  in  an  equivalent  Hugoniot  elastic  limit  for 
plane  strain,  cthel,ps,  given  by  Equation  2.15.  Due  to  the  limitations  of  strain  in 
the  uniaxial  strain  case  and  the  plane  strain  case,  it  was  expected  that  yield  stress, 
Ya,  would  be  less  than  both  <7hel,ps  and  &hel-  Also,  that  <Jhel,ps  would  be  less 
than  ctHEl ■  The  assumption  was  validated  when  the  equations  derived  and  solved  for 
VascoMax  300  in  Table  2.1.  Equation  2.16  was  given  in  this  chapter  as  a  method  to 
solve  for  the  speed  of  an  elastic  wave  through  a  solid  material  under  plane  strain  con¬ 
ditions.  Appendix  A  includes  the  full  derivation  of  the  equations  used  to  determine 
the  equivalent  Hugoniot  elastic  limit  and  plane  strain  elastic  wave  speed. 

Most  of  the  previous  work  used  a  finite  element  code  to  simulate  the  problem. 
Meador  used  CTH,  a  hydrocode  discussed  in  Section  2.8,  to  model  the  wear  of  HHSTT 
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slippers.  This  research  uses  the  same  code  used  by  Meador,  with  some  modifications 
mentioned  in  Chapter  III.  This  chapter  discussed  in  detail  the  benefits  of  using  a 
hydrocode  as  opposed  to  a  finite  element  code.  A  primary  benefit  includes  the  use  of 
an  Eulerian-Lagrangian  coordinate  system  to  avoid  large  mesh  distortions,  and  allow 
for  high  velocity  impact  scenarios. 

Hale  [20]  conducted  metallurgical  studies  on  both  used  and  unused  test  slippers 
from  the  HHSTT.  The  study  suggested  that  mechanical  wear  results  from  plastic 
deformation.  This,  along  with  the  fact  that  the  micro-level  simulation  technique, 
first  proposed  by  Chmicl,  is  a  time-dependent  process  requires  the  use  of  a  viscoplas¬ 
tic  constitutive  model.  The  Johnson-Cook  model,  Equation  2.7  was  chosen  for  this 
research,  because  it  includes  considerations  of  large  strains,  high  strain  rates,  and  ele¬ 
vated  temperatures.  Furthermore,  the  equation  was  intended  for  use  in  computational 
software. 

Although  wear  involves  the  removal  of  material,  developing  a  model  to  remove 
material  during  simulation  would  be  complicated.  Therefore,  a  qualitative  measure, 
based  on  material  damage,  has  been  adopted.  The  Johnson-Cook  constitutive  model 
was  used  to  develop  failure  criteria,  discussed  in  Section  2.9.  The  failure  criterion 
is  used  to  evaluate  material  as  damaged  or  undamaged.  The  amount  of  damage 
material  for  each  simulation  is  computed  and  divided  by  the  distance  slid  to  give  a 
wear  rate.  Previous  work  by  Hale  and  Meador  have  made  use  of  models  developed  by 
Archard  [4]  to  relate  the  two-dimensional  plane  strain  single  asperity  collision  event 
to  a  three  dimensional  wear  scenario. 
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III.  Numerical  Model 


Developing  a  numerical  model  to  simulate  entire  HHSTT  missions  is  an  imprac¬ 
tical  approach  to  the  wear  phenomenon,  in  terms  of  run-time  and  simulation  cost. 
However,  previous  research  by  Chrniel  [14],  Burton  [10],  Hale  [20],  Meador  [28],  and 
Lodygowski  [27]  has  shown  results  using  a  simplified  plane-strain  model.  The  models 
in  previous  research  have  made  use  of  the  DADS  data  provided  by  the  HHSTT.  This 
chapter  will  discuss  the  DADS  system  and  how  the  recorded  data  is  used  to  character¬ 
ize  the  slipper-rail  sliding  event.  This  chapter  will  also  present  the  hydrocode  model 
used  for  this  research,  discussing  the  input  parameters  including  initial  velocity,  the 
viscoplasticity  model,  and  equation  of  state.  The  method  used  to  calculate  plane 
strain  mechanical  wear  rates,  and  total  mechanical  wear  will  also  be  presented. 

3.1  DADS  Data 

Properly  characterizing  the  slipper  dynamics  as  it  slides  along  the  rail  is  neces¬ 
sary  to  create  an  accurate  model.  The  HHSTT  provides  data  using  a  program  called 
Dynamic  Analysis  and  Design  System  (DADS).  DADS  is  a  commercial-off-the-shelf 
software  developed  by  Computer  Aided  Design  Software,  Inc.  The  HHSTT  uses  the 
software  to  simulate  a  rocket  sled  run  and  predict  vertical  forces  of  each  slipper,  ver¬ 
tical  velocity  of  the  sled,  and  horizontal  velocities  of  the  slippers,  as  a  function  of 
time. 

A  model  has  been  developed  to  represent  the  HHSTT  sled  and  rail  as  a  com¬ 
plicated  system  of  masses,  springs,  and  dampers,  while  the  sled  forward  velocity  and 
rail  undulations  are  supplied  as  inputs  to  the  system  [21],  The  simulations  have  been 
validated  using  accelerometers  attached  to  test  sleds  [22],  Given  the  complexity  of 
the  model,  and  the  validation  by  HHSTT  engineers,  the  DADS  model  is  assumed  to 
be  valid  within  the  context  of  this  research  [28]. 

The  geometry  of  the  sled  influences  the  dynamics  of  the  sled  and  slipper  as  it 
slides  along  the  steel  rail.  According  to  the  HHSTT  Design  Manual  [1],  the  nominal 
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slipper  gap  is  0.125  inches.  This  corresponds  to  the  nominal  max  clearance  between 
the  slipper  and  rail.  The  slipper  gap  allows  the  slipper  to  bounce  along  the  rail 
during  the  test  missions.  Due  to  this  gap  and  bouncing  effect,  the  slipper  is  not  in 
total  contact  with  the  rail  for  the  entire  run.  The  bouncing  effect  is  included  in  the 
calulation  of  total  mechanical  wear  in  Section  3.4. 


3.2  Plane  Strain  Simulation  Using  a  Hydrocode 

The  plane  strain  simulations  for  this  thesis  used  CTH,  a  hydrocode  discussed  in 
Section  2.8.  The  model  simulates  the  collision  of  VascoMax  300  with  a  6  /mi  radius 
hemispherical  surface  asperity  made  of  AISI  1080  steel,  as  shown  in  Figure  3.1.  In  this 
figure,  the  slipper  (yellow)  moves  to  the  right  at  a  given  input  velocity  and  collides 
with  the  green  asperity  and  rail.  Several  inputs  must  be  defined  to  run  the  two- 
dimensional  analysis.  These  inputs  include  sliding  velocities,  simulation  time,  mesh 
and  domain  sizing,  boundary  conditions,  and  geometry.  This  research  made  use  of  an 
existing  CTH  code  developed  by  Meador  [28]  with  some  modifications. 

The  sliding  distance  was  chosen  to  be  110%  of  the  6  /mi  radius  to  allow  the 
leading  edge  of  the  slipper  to  go  past  the  maximum  height  of  the  asperity.  The 
simulation  time  is  found  as  a  function  of  the  input  sliding  velocity,  vsude,  and  asperity 
radius,  ra.  Equation  3.1  was  used  to  determine  the  simulation  time  for  each  run.  The 
simulation  times  are  shown  in  Table  3.1. 


,  _  (l-l)(r«) 

t'sim 

V  slide 

The  size  of  the  domain  and  slipper  were  selected  to  reduce  boundary  effects  and 
pressure  wave  interactions  along  the  edges.  Since  the  simulation  time  is  known  for 
each  case,  and  the  velocity  of  the  elastic  pressure  wave  and  plastic  pressure  waves  are 
given  by  Equations  (2.13)  and  (2.14),  the  distance  traveled  by  a  pressure  wave  can  be 
found.  Meador  determined  that  a  domain  size  of  850  /mi  by  850  /mi,  and  a  slipper  size 
of  approximately  700  /mi  by  125  /mi  were  sufficient.  However,  for  his  research,  Meador 


(3.1) 


43 


0.0215 


0.021 


-v  0.0205 

rH 

e 

CJ 

^  0.02 


0.0195 


0.019 


Materials  at  6.60e-ll  seconds 

i — 1 — 1 — ' — ■ — i — ' — ■ — ' — ’ — i— 


0.069  0.07  0.071 


X  (cm) 

1080  Steel  Rail 
VascoMax  300  Slipper 

Figure  3.1:  Slipper- Asperity  Interface  of  Current  Model 
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Table  3.1:  Simulation  Time  for  Given  Velocity 

Based  on  Equation  3.1 


Horizontal  Velocity  (m/s) 

Simulation  Time  (s) 

100 

6.60  x  10"8 

200 

3.30  x  10"8 

300 

2.20  x  10"8 

400 

1.65  x  10"8 

500 

1.32  x  10"8 

600 

1.10  x  10"8 

700 

9.43  x  10"9 

800 

8.25  x  10"9 

900 

7.33  x  10"9 

1,000 

6.60  x  10"9 

1,100 

6.00  x  HT9 

1,200 

5.50  x  10"9 

1,300 

5.08  x  10"9 

1,400 

4.71  x  10"9 

1,500 

4.40  x  10"9 

considered  velocities  ranging  from  750  m/s  to  3,000  m/s.  This  research  considers  a 
velocity  range  from  200  m/s  to  1,500  m/s,  corresponding  with  the  velocity  profile 
of  the  third  sled  from  the  January  2008  test  mission  (1.12).  Decreasing  the  sliding 
velocity  increases  the  simulation  time,  Equation  3.1,  which  increases  the  distance 
traveled  by  the  stress  waves.  This  means  at  velocities  less  than  750  m/s,  the  stress 
waves  may  reach  the  boundaries  of  the  slipper  or  domain. 

To  reduce  the  wave  interactions  along  the  boundary  of  the  domain,  the  boundary 
conditions  were  modified.  A  combination  of  Type  1  and  Type  2  boundary  conditions 
were  used  to  simulate  a  semi-infinite  boundary  (Type  1  and  2)  and  allow  material  to 
flow  out  of  the  domain  (Type  2  only).  Figure  3.2  shows  the  domain  created  with  the 
materials  defined  and  boundary  conditions  selected. 

Additional  consideration  was  given  to  the  leading  edge  of  the  slipper,  specifically 
at  the  point  of  contact  with  the  surface  asperity.  The  concern  was  that  the  corner 
would  result  in  a  singularity.  Previous  research  by  Cameron  [13],  Hale  [20],  and 
Meador  [28]  alleviated  this  issue  by  replacing  the  corner  with  a  2  /mi  radius  fillet. 
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Materials  at  3.30e-10  seconds 
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Figure  3.2:  Entire  Domain  of  Current  Model 
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This  research  uses  the  filleted  edge  to  stay  consistent  with  the  previous  work.  The 
mesh  was  selected  to  create  a  grid  of  cells  with  size  1  //rn  by  1  /irn  throughout  the 
entire  domain.  Figure  3.3  shows  the  mesh  used  for  the  simulations. 
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Materials  at  1.80e-10  seconds 
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Figure  3.3:  Eulerian  Mesh  Applied  to  Current  Model 


The  tracer  points,  discussed  in  Section  2.8.4  are  attached  to  the  material  through¬ 
out  the  entire  simulation.  This  method  of  data  collection  allows  each  cell  to  be  affected 
by  the  previous  iteration,  meaning  if  a  cell  has  failed  at  one  time  step  it  will  remained 
failed  for  the  rest  of  the  simulation.  The  tracer  points  are  initially  placed  at  the  center 
of  each  cell,  so  the  area  allocated  to  each  tracer  is  the  same  as  the  cell  area.  MATLAB 
was  used  to  create  a  post-processing  code  to  compute  the  wear  rates  from  each  CTH 
run.  The  post-processing  code  is  provided  with  explanation  in  Appendix  C. 


3.2.1  Material  Interface  Conditions.  The  model  developed  in  this  thesis 
allows  one  material  to  slide  along  another  surface.  Therefore,  friction  along  the  inter¬ 
face  must  be  considered.  However,  CTH  does  not  allow  an  input  coefficient  of  friction. 
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There  are  two  interface  conditions  that  can  be  defined  in  the  model  to  represent  fric¬ 
tion  in  the  surface.  The  first  is  called  a  “slide  line”  condition  which  sets  the  shear 
stress  along  the  interface  to  zero.  This  condition  is  used  to  simulate  a  frictionless 
surface.  The  second  condition  is  called  “no  slide”.  This  condition  will  not  allow  the 
material  to  move  until  a  user  defined  pressure  is  reached.  This  user  defined  pressure  is 
the  fracture  pressure  defined  in  the  fracture  input  set  of  the  CTH  code  in  Appendix  B. 
The  CTH  code  was  run  with  these  two  conditions  representing  two  extremes  in  terms 
of  friction;  one  being  a  frictionless  surface,  and  the  other  representing  a  semi-infinite 
coefficient  of  friction. 

3.2.2  Input  Velocity.  Since  the  simulation  is  evaluating  a  two-dimensional 
scenario,  velocity  can  be  defined  in  two  directions:  horizontal  and  vertical.  The 
horizontal  component  of  velocity  is  determined  by  the  sliding  velocity.  Values  are 
chosen  from  200  m/s  to  1,500  m/s  to  represent  the  increasing  velocity  of  the  slipper 
along  the  rail.  The  vertical  velocity  component  can  be  determined  using  the  DADS 
data.  The  vertical  velocity  is  positive  going  up  and  negative  going  down.  Therefore, 
the  negative  represents  the  slipper  moving  toward  the  rail.  Figure  3.4  plots  the 
vertical  velocity  of  the  third  stage  aft  right  slipper  from  the  2008  mission  against 
the  horizontal  velocity  of  the  sled.  The  figure  shows  that  the  maximum  vertical 
velocity  of  the  slipper  into  the  rail  is  approximately  -3.45  m/s.  The  horizontal  velocity 
component  is  always  much  larger  than  the  vertical  component.  This  means  that 
vertical  velocity  has  very  little  effect  on  the  velocity  vector.  All  simulations  have  a 
constant  vertical  velocity  of  -0.5  m/s,  which  helps  to  keep  the  slipper  in  contact  with 
the  rail  throughout  the  simulation.  The  idea  of  creating  a  more  realistic  simulation  of 
the  HHSTT  environment  is  continued  with  the  addition  of  a  dead  load  to  represent 
the  effective  mass  of  the  slipper  as  it  slides  along  the  rail  at  increasing  velocity.  The 
addition  of  the  dead  load  is  discussed  in  detail  in  Section  4.1. 

3.2.3  Viscoplasticity  Model.  The  Johnson-Cook  constitutive  equation,  dis¬ 
cussed  in  Section  2.4,  is  used  as  the  viscoplasticity  model  for  this  research.  The  model 
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Aft  Rear  Slipper  Vertical  Velocity  vs.  Sled  Forward  Velocity  Maximum  Donward  Velocity  =  -3.45  m/s 
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Figure  3.4:  Slipper  Vertical  Velocity  from  DADS 


was  developed  to  handle  large  strains,  high  strain  rates,  and  high  temperatures.  The 
material  constants  need  to  be  defined  for  the  VascoMax  300  slipper  and  AISI  1080 
steel  rail  in  order  to  use  the  Johnson-Cook  model.  Cinnamon  [15-17]  determined 
these  constants,  shown  in  Table  3.2,  using  flyer  plate  experiments.  It  is  also  necessary 
to  define  the  initial  temperature  of  the  materials  because  the  Johnson-Cook  model 
carries  the  homologous  temperature  term,  Equation  2.8. 


Table  3.2:  Johnson-Cook  Coefficients  for  VascoMax  300  and  AISI  1080  Steel 


Coefficient 

VascoMax  300 

AISI  1080  Steel 

A  (GPa) 

2.1 

0.7 

B  (GPa) 

0.124 

3.6 

C  (Unitless) 

0.03 

0.17 

m  (Unitless) 

0.8 

0.25 

n  (Unitless) 

0.3737 

0.6 

Tmelt  (K) 

1,685 

1,630 

3.2.4  Equation  of  State.  Hydrocodes  make  use  of  an  equation  of  state  to 
relate  internal  energy  and  pressure  of  a  material  to  the  density  and  temperature.  It 
serves  as  an  additional  equation  to  relate  the  conservation  equations  to  the  consti¬ 
tutive  equation.  The  equation  of  state  is  also  useful  when  a  shock  is  present  within 
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the  material.  The  shock  creates  discontinuities  and  an  EOS  can  be  used  to  solve  for 


material  properties. 

The  Mie-Griineisen  EOS,  presented  in  Section  2. 8. 2.1,  was  initially  considered 
for  use  in  this  research.  This  EOS  model  is  typically  used  for  high  velocities  ranging 
from  500  m/s  to  2,000  m/s  [37].  Therefore,  some  issues  arose  when  modeling  at  the 
lower  sliding  velocities.  These  errors  included  numerical  inconsistencies  in  calculated 
mechanical  wear  rates  and  pressure  wave  propagation.  These  low  velocity  issues  and 
a  modified  approach  are  discussed  in  greater  detail  in  Section  4.5. 

Due  to  the  inconsistent  results  using  the  Mie-Griineisen  equation  of  state  at  low 
velocities,  the  Sesame  EOS  is  used  for  this  research.  The  Sesame  EOS  is  a  tabular  set 
of  experimental  data.  The  experiments  are  typically  a  high  velocity  impact  scenario 
under  uniaxial  strain  conditions.  VascoMax  300  is  defined  in  the  Sesame  tables  and  is 
used  to  represent  the  slipper  for  this  research.  However,  AISI  1080  steel  is  not  defined 
in  the  Sesame  tables,  so  iron  is  used  to  represent  the  rail.  The  two  materials  have 
similar  properties  as  shown  in  Table  3.3. 


Table  3.3:  Iron  and  AISI  1080  Steel  Properties 


Property 

Iron 

AISI  1080  Steel 

Density  (g/cm3) 

7.28 

7.85 

Yield  Stress  (MPa) 

50 

585 

Elastic  Modulus  (GPa) 

200 

205 

Melt  Temperature  (K) 

1,181 

1,630 

Poisson’s  Ratio 

0.28 

0.25 

3.3  Mechanical  Wear  Rate  Calculation 

The  post-processing  code  uses  CTH  output  to  determine  the  wear  rate  per  unit 
width,  Wuw,  given  by  Equation  3.2,  where  Ad  is  the  damage  area  computed  from  the 
plane  strain  simulation  based  on  the  failure  criteria  used,  and  the  distance  slid  is  the 
product  of  sliding  velocity,  vsude,  and  simulation  time,  tsim. 
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Wuw  =  Ad  (3.2) 

^slide^sim 

The  plane  strain  simulation,  along  with  this  equation  gives  wear  rates  in  units 
of  area  of  damaged,  or  worn,  material  per  distance  slid.  Since  wear  is  defined  as  the 
volume  of  material  worn  per  distance  slid,  it  is  best  represented  as  a  three-dimensional 
problem.  As  such,  a  conversion  factor  must  be  established  to  represent  a  three- 
dimensional  hemispherical  surface  asperity  using  the  two-dimensional  semi-circular 
asperity  in  plane  strain. 

3.3.1  Semi- spherical  Coefficient.  Hale  [20]  determined  the  semi- spherical 
coefficient  by  running  plane  strain  simulations  with  2  /nn,  4  //rn,  and  6  /nn  asperities 
and  integrating  across  the  width  to  determine  a  volume  of  damaged  material  per 
distance  slid.  Figure  3.5  shows  how  the  plane  strain  models  are  related  to  the  three- 
dimensional  analysis.  The  red  areas  in  Figure  3.5(b)  represent  the  damaged  material 
due  to  the  collision  with  the  various  asperities.  For  a  given  sliding  velocity,  the  area 
of  damaged  material  increases  as  the  size  of  the  asperity  increases. 

The  2  /mi  and  4  /nn  asperity  collisions  are  related  to  the  6  /mi  asperity  by 
assuming  an  off-center  collision.  Since  the  analysis  is  under  plane  strain  conditions 
and  the  z-axis  is  eliminated,  the  actual  height  of  the  asperity  does  not  affect  the 
simulation.  Equation  3.3  is  used  to  determine  the  location  of  the  2  /nn  and  4  /nn 
asperity  along  the  z-axis  in  the  three-dimensional  hemispherical  asperity,  where  r  is 
the  6  /nn  asperity  radius. 


z  =  \Jr2  -  y 2  (3.3) 

This  places  the  2  /nn  and  4  /nn  asperities  at  z-locations  of  4.47  /nn  and  5.66 
/nil,  respectively.  These  locations  are  illustrated  by  the  dashed  line  in  Figure  3.5(a). 
The  single  asperity  wear  rate,  Wsa,  is  determined  by  Equation  3.4,  where  the  integral 
is  multiplied  by  two  to  represent  the  symmetrical  asperity. 
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(a)  Hemispherical  Asperity 


6  pm  4  pm  2  pm 


(b)  Plane  Strain  Analyses  Illustration 


Figure  3.5:  Plane  Strain  Representation  of  a  Semi-Spherical 
Surface  Asperity  [20] 
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Wuw(z)dz 


(3.4) 


WH, 


Equation  3.4  gives  the  plane  strain  wear  rate  for  a  collision  with  a  single  semi- 
spherical  asperity.  Hale  solved  this  equation  for  a  range  of  sliding  velocities,  and  then 
divided  the  single  asperity  wear  rates  by  their  respective  plane  strain  wear  rates  to  get 
an  average  effective  width,  we//,  of  8.29  /irn.  Equation  3.5  uses  this  average  effective 
width  to  calculate  the  single  asperity  wear  rates,  as  opposed  to  using  the  integral  in 
Equation  3.4. 


bEs  a  —  '^e//bEl[UI  (3-5) 

3.3.2  Archard  Scaling  Factor.  The  plane  strain  models  developed  by  Hale 
[20]  and  Meador  [28]  make  use  of  a  scaling  factor  to  account  for  collisions  with  multiple 
asperities  as  the  slipper  sides  along  the  rail.  This  equation  is  derived  by  relating  wear 
rates  to  Archard’s  wear  model  at  low  velocities  [2-4],  Equation  3.6  is  used  to  relate 
Archard’s  wear  rate,  W a ,  to  the  single  asperity  wear  rate. 

wA=kjf  =  NWsa  (3.6) 

where  Ia  is  Archard’s  wear  coefficient,  F  is  the  applied  load,  FI  is  the  material  hard¬ 
ness,  and  N  is  the  scaling  factor.  The  applied  load  in  Archard’s  equation  relates  to 
the  force  applied  by  the  pin  in  a  pin-on-disk  experiment.  The  scaling  factor  is  found 
by  Equation  3.7.  Hale  [20]  solved  for  N  =  11.77,  at  a  sliding  velocity  of  10  m/s,  with 
Ia  equal  to  4.45  x  10-5,  and  with  F  given  from  the  DADS  data. 


N 


kAF 

WxnH 


(3.7) 
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3-4  Total  Mechanical  Wear  Calculation 

Calculating  total  wear  of  the  HHSTT  slippers  allows  for  a  comparison  to  be 
made  between  experimental  data  and  simulation  results.  The  experimental  data  for 
this  research  is  a  used  slipper  retrieved  from  the  January  2008  test  mission  at  Hollo¬ 
man  AFB.  To  evaluate  total  wear,  the  single  asperity  wear  rates  are  integrated  with 
respect  to  sliding  distance.  The  DADS  data  is  used  to  determine  at  what  distance 
along  the  track  the  slipper  reaches  a  certain  velocity.  The  wear  rates  are  then  plot¬ 
ted  as  a  function  of  sliding  distance.  Before  the  values  are  integrated,  an  additional 
scaling  factor  must  be  included. 

This  scaling  factor  represents  the  amount  of  time  the  slipper  and  rail  are  in 
contact.  The  CTH  simulations  assume  the  slipper  and  rail  are  in  contact  throughout 
the  entire  simulation.  The  percentage  of  contact,  dpc,  is  set  to  0.3.  That  assumes 
that  the  slipper  and  rail  are  in  contact  for  30%  of  the  test  run.  This  assumption,  first 
used  by  Hale  [20],  comes  from  the  comparison  of  test  data  between  the  January  2008 
mission  and  a  simulation  designated  by  the  HHSTT  as  80X-A1.  DADS  data  for  the 
80X-A1  simulation  was  supplied  to  Cameron  [13]  for  his  research  in  2007. 

Equation  3.8  calculates  total  mechanical  wear,  Wtotal,  where  dmax  is  the  to¬ 
tal  sliding  distance.  This  equation  includes  the  single  asperity  wear  rate,  Wsa  from 
Equation  3.5,  the  Archard  scaling  factor,  N ,  from  Equation  3.7,  and  the  percentage 
of  contact  coefficient,  dpc. 


["dmax 

Wtotal  =  Ndpc  /  Wsa(s)ds  (3.8) 

Jo 

3. 5  Summary  of  Numerical  Modeling 

This  chapter  discussed  the  dynamic  data,  DADS,  used  to  characterize  position, 
velocity,  and  forces  of  the  rocket  sled  system.  This  dynamic  data,  along  with  the 
theoretical  background  presented  in  Chapter  II,  was  used  to  create  a  hydrocode  model 
capable  of  estimating  mechanical  wear  of  HHSTT  slippers.  Equations  were  given  in 
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this  chapter  to  represent  multiple  semi-spherical  asperity  collisions  with  a  single  semi¬ 
circular  asperity  in  plane  strain.  The  model  developed  in  this  thesis  evaluates  pressure 
and  internal  energy  of  a  material  due  to  this  collision.  The  next  chapter  will  discuss 
the  results  of  the  simulation,  calculated  mechanical  wear  rates,  total  mechanical  wear, 
and  compare  these  to  the  experimental  data. 
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IV.  Results  and  Discussion 


This  chapter  will  present  the  results  of  the  numerical  model  discussed  in  Chap¬ 
ter  III.  First,  the  equivalent  plane  strain  Hugoniot  clastic  limit,  discussed  in  Sec¬ 
tion  2.5.2,  will  be  compared  against  CTH  simulations.  The  calculated  wear  rates  will 
be  presented  for  various  failure  criteria  and  interface  conditions.  Finally,  the  results 
of  the  total  wear  calculation  will  be  presented  and  compared  against  experimental 
wear  from  the  HHSTT  January  2008  test  mission. 

4-1  Dead  Load 

In  order  to  develop  a  model  that  represents  the  HHSTT  environment,  a  variable 
vertical  force  was  considered.  Since  the  concept  of  a  hydrocode  is  to  include  kinetic 
energy  as  a  basic  function,  it  was  felt  that  at  least  an  investigation  of  the  vertical 
movement  should  be  included.  This  is  considered  by  characterizing  a  vertical  force. 
This  variable  force  represents  the  vertical  force  of  the  slipper  as  it  bounces  along  the 
rail.  Figure  4.1  is  a  plot  of  the  vertical  force  of  the  aft  right  slipper  from  the  third 
sled  during  the  January  2008  test  from  the  DADS  data.  When  the  force  is  zero,  the 
slipper  is  not  in  contact  with  the  rail. 

CTH  does  not  allow  a  force  input.  However,  the  force  can  be  represented  with 
an  appropriate  dead  load  fixed  to  the  top  of  the  slipper.  DADS  does  not  provide  data 
representing  vertical  acceleration  of  the  slipper  as  it  travels  down  the  rail.  Therefore, 
a  modified  force  equation,  Equation  4.1,  was  developed  to  represent  the  effective  dead 
load  as  a  function  of  vertical  velocity  and  vertical  force  using  Newton’s  second  law, 
F=ma. 


m  = 


F  At 

■ L  avg 

Av 


(4.1) 


where  Av  is  the  maximum  change  in  vertical  velocity  within  a  window  enclosing  the 
desired  horizontal  velocity,  At  is  the  size  of  the  window,  and  Favg  is  the  average 
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Aft  Right  Slipper  Contact  Forces 


Figure  4.1:  HHSTT  Third  Stage  Vertical  Force,  January  2008 
Test  Mission 


slipper  force  during  the  window.  The  velocity  profile  of  the  third  stage  from  the 
HHSTT  January  2008  test  mission,  Figure  4.2,  is  used  to  determine  the  time  at 
which  a  desired  velocity  is  reached.  Figures  4.3  and  4.4  show  the  method  of  applying 
windows  surrounding  a  time  representing  a  target  sled  velocity  of  500  m/s  enclosing 
the  vertical  force  and  vertical  velocity. 

The  size  of  the  window,  At,  is  the  time  between  two  peaks  enclosing  the  desired 
velocity.  These  two  peaks  are  shown  in  Figure  4.3  with  red  circles.  The  black  circle 
represents  the  time  at  which  the  sled  reaches  the  desired  velocity.  The  average  force, 
Favg,  is  the  average  of  the  two  peaks  in  Figure  4.3.  The  same  window  applied  to 
the  vertical  force  plot  is  applied  to  the  vertical  velocity.  A  maximum  and  minimum 
velocity  is  found  in  this  window  and  shown  on  Figure  4.4  with  two  red  circles.  The 
maximum  change  in  velocity,  Av  is  the  change  in  velocity  between  these  two  points. 

The  dead  load  is  calculated  with  respect  to  sliding  velocity  and  plotted  in  Fig¬ 
ure  4.5.  It  is  important  to  note  that  there  are  many  sudden  changes  in  vertical  velocity 
and  vertical  force  of  the  slipper  as  it  travels  along  the  rail.  The  scattered  data  results 
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Sled  C.G.  Velocity 


Figure  4.2:  HHSTT  Third  Stage  Velocity  Profile,  January 

2008  Test  Mission 


x  10“ 


time  (s) 

Figure  4.3:  Windowed  Vertical  Force  Data  at  500  m/s 
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Figure  4.4:  Windowed  Vertical  Velocity  Data  at  500  m/s 

in  a  scattered  function  of  dead  load  with  respect  to  sliding  velocity.  Also,  there  is 
an  outlying  maximum  dead  load  at  1,000  m/s.  This  is  probably  due  to  the  increased 
acceleration  down  the  track  due  to  the  bring  of  the  third  rocket  sled. 

The  model  developed  in  this  thesis  is  a  local  submodel  of  the  slipper  colliding 
with  a  surface  asperity.  As  such,  the  entire  slipper  is  not  represented  in  the  model. 
Therefore,  a  mass  fraction,  5m,  was  found  that  related  the  mass  of  the  actual  slipper 
to  the  mass  of  the  slipper  in  the  CTH  model.  This  mass  fraction  is  applied  to  the 
calculated  dead  load  when  added  to  the  CTH  input  deck.  Data  from  Holloman  AFB 
states  that  the  weight  of  the  aft  right  slipper  is  19  pounds.  This  is  converted  to  grams, 
assuming  gravity,  g,  =  9.81  m/s2  using  Equation  4.2.  This  equation  gives  a  slipper 
mass  of  approximately  8,615  grams.  The  mass  of  the  CTH  slipper  is  given  by  the 
product  of  its  density,  8.091  g/ern3  and  volume,  8.63  xlO-4  cm3,  to  be  approximately 
7.02  xl0~3  grams.  Therefore,  the  mass  fraction,  5m,  is  8.15  x  10~'. 


Mass 


grams 


.  4.448V  kg  ■  m  1  1000g 

ei9  hbf  jy  _  s2  g.3im/s2  Ykg 


(4.2) 
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Figure  4.5:  Change  in  Dead  Load  with  Respect  to  Velocity 

To  add  the  calculated  dead  load  to  the  CTH  model,  a  block  of  platinum  was 
fixed  to  the  top  of  the  VascoMax  300  slipper.  Platinum  was  selected,  because  it  is 
the  most  dense  material  in  CTH.  Since  the  dead  load  sits  atop  the  entire  slipper,  the 
dead  load  dimension  along  the  X-axis  is  known.  Because  a  plane  strain  model  was 
developed,  the  thickness  is  also  known.  Since  the  mass  of  any  material  is  given  by 
the  product  of  its  density  and  volume,  the  height  of  the  platinum  dead  load  is  found 
using  Equation  4.3,  where  h  is  the  dead  load  height,  m  is  the  dead  load  mass,  p  is 
the  dead  load  density,  w  is  the  width  of  the  dead  load  along  the  X-axis,  and  t  is  the 
thickness  of  the  dead  load  in  the  z-axis. 


mSM 

h  = - 

pwt 


(4.3) 


The  investigation  of  the  dead  load  effect  was  carried  out  by  running  two  sim¬ 
ulations  with  the  dead  load  at  800  m/s  and  1,200  m/s  with  a  “no  slide”  boundary 
condition.  The  results  of  these  simulations,  shown  in  Table  4.1,  are  identical  to  the 
simulations  without  the  dead  load.  This  is  probably  due  to  the  fact  that  the  con¬ 
tact  area  between  the  slipper  and  asperity  is  so  small  that  the  added  mass  does  not 
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have  a  significant  effect.  The  dead  load  was  dropped  from  the  simulation,  and  is  not 
considered  in  any  more  simulations. 


Table  4.1:  Dead  Load  Wear  Rates 


Horizontal  Sliding  Velocity  (m/s) 

800 

1,200 

Dead  Load  Height  (cm) 

4.30xl0~2 

3.04xl0-2 

Dead  Load  Strain  at  Max  Stress  (mm3/mm) 

2.66xl0-4 

2.86X10"4 

Dead  Load  Critical  Von  Mises  Stress  (mm3/mm) 

3.54xl0-4 

3.73xl0~4 

No  Dead  Load  Strain  at  Max  Stress  (mm3 /mm) 

2.66xl0~4 

2.86X10"4 

No  Dead  Load  Critical  Von  Mises  Stress  (mm3/mm) 

3.54xl0“4 

3.73xl0-4 

4-2  Failure  Criteria  Selection 

Three  failure  Criteria  were  presented  in  Section  2.9,  plastic  strain  at  max  stress, 
critical  von  Mises  stress,  and  the  Johnson-Cook  fracture  model.  The  plastic  strain  at 
max  stress  criteria,  developed  by  Meador  [28],  has  provided  reliable  results  in  previous 
research  and  is  considered  a  valid  failure  criteria.  The  critical  von  Mises  stress  criteria 
is  a  modified  approach  used  by  Hale  [20]  and  Meador.  The  benefit  of  this  criteria  is  in 
the  fact  that  the  von  Mises  stress  can  be  calculated  from  the  CTH  simulation.  This 
removes  the  need  for  the  curve  fit,  Equation  2.22.  The  number  of  cells  exceeding  the 
critical  von  Mises  stress  value  is  used  to  determine  the  amount  of  material  damage 
for  this  criteria.  A  critical  stress  value  of  3,000  MPa  was  selected  based  upon  Hale’s 
strain  rate  analysis  [20]. 

Preliminary  evaluation  of  the  Johnson-Cook  fracture  model  suggest  that  it  does 
not  work  well  with  the  plane  strain  collision  studied  in  this  research.  This  method 
produced  zero  wear  in  most  cases.  Recall  that  the  Johnson-Cook  fracture  model  is 
dependent  upon  an  integral  given  by  Equation  4.4.  As  damage  accumulates  in  the 
material,  the  integral  goes  to  1.  When  the  integral  is  equal  to  1,  the  material  is  said  to 
have  failed.  For  most  cases,  the  model  developed  for  this  research  does  not  allow  the 
integral  to  reach  the  critical  value  of  1.  However,  some  cells  do  reach  the  critical  value 
resulting  in  some  wear.  This  number  tends  to  be  two  to  three  orders  of  magnitude 
less  than  the  other  models. 
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dep 

£Pf(p,  Y,  T,  e) 


(4.4) 


D  = 


One  seemingly  obvious  reason  for  the  poor  results  of  the  model  is  the  fact  that 
the  coefficients  of  the  Johnson-Cook  fracture  model  for  VascoMax  300  and  AISI  1080 
steel  are  not  defined  in  CTH.  Several  steels  are  defined  in  CTH.  For  the  purposes 
of  this  research.  Iron  was  used  to  represent  the  AISI  1080  steel  rail,  and  AISI  4340 
steel  was  used  to  represent  the  VascoMax  300  slipper.  The  Johnson-Cook  fracture 
coefficients  used  are  given  in  Table  4.2.  Previous  work  by  Lee  [25]  has  shown  the  utility 
of  the  Johnson-Cook  fracture  model  in  CTH.  If  the  coefficients  for  the  materials  were 
defined  in  CTH,  the  model  may  produce  satisfactory  results. 


Table  4.2:  Johnson-Cook  Fracture  Coefficients  for 

Iron  and  AISI  4340  Steel  Defined  in  CTH  [24] 


Coefficient 

Iron 

AISI  4340  Steel 

D\ 

-2.2 

-0.8 

d2 

5.43 

2.1 

d3 

-0.47 

-0.5 

d4 

0.016 

0.002 

d5 

0.63 

0.61 

Tmelt  (eV) 

0.1581885 

0.1566566 

4-3  Validation  of  Plane  Strain  Hugoniot  Limit 

Table  2.1  provides  the  predicted  equivalent  plane  strain  Hugoniot  elastic  limit 
for  VascoMax  300  as  2.8664  GPa.  This  value  can  be  validated  by  checking  for  an 
equivalent  Hugoniot  limit  from  the  CTH  simulation.  This  is  done  by  plotting  the 
evolution  of  pressure  through  time  at  a  point  in  the  VascoMax  300  slipper.  Cinnamon 
[15]  did  this  under  uniaxial  strain  conditions  using  a  flyer  plate  test.  The  experiment 
bred  VascoMax  300  projectiles  at  a  target  at  high  velocites.  Stress  was  measured 
using  a  stress  gauge  attached  to  the  projectile  approximately  2  mm  from  the  leading 
edge.  Since  the  experiment  represented  uniaxial  strain,  the  pressure  was  set  to  the 
measured  stress.  The  pressure  was  plotted  against  time  to  check  for  the  Hugoniot 
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clastic  limit.  Cinnamon  also  created  a  CTH  simulation  of  the  impact  event.  To 
keep  in  line  with  the  experiment,  the  pressure  from  Cinnamon’s  CTH  simulation  was 
recorded  2  mm  from  the  interface  of  the  projectile  and  target.  The  dimensions  of 
the  slipper-rail  simulation  created  for  this  thesis  are  much  smaller  than  Cinnamon’s 
model.  Therefore,  the  pressure  is  recorded  from  the  CTH  simulation  at  2  /jrn  vertically 
from  the  interface  of  the  slipper  and  rail,  and  2  /im  horizontally  from  the  interface  of 
the  slipper  and  asperity.  The  black  dot  in  Figure  4.6  shows  the  point  in  the  model 
where  the  pressure  data  is  recorded  to  locate  the  equivalent  plane  strain  Hugoniot 
clastic  limit.  Figure  4.7  shows  the  increase  of  pressure,  in  a  slipper  sliding  at  1,000 
m/s,  to  a  value  in  which  the  pressure  drops  and  increases  to  a  max  value.  The  value 
at  which  the  pressure  stops  increasing  is  referred  to  as  the  equivalent  plane  strain 
Hugoniot  clastic  limit.  When  loading  exceeds  this  value,  deformation  is  no  longer 
purely  clastic.  Multiple  plastic  waves  are  produced  when  the  loading  exceeds  the 
limit,  which  explains  the  cyclic-  behavior  at  max  pressure  after  the  equivalent  plane 
strain  Hugoniot  limit  is  reached. 

It  should  be  noted  that  the  equivalent  plane  strain  Hugoniot  limit  from  Fig¬ 
ure  4.7  does  not  equal  the  calculated  value  (2.8664  GPa)  exactly.  However,  the  value 
from  the  simulation  is  close  to  the  predicted  value.  Figures  4.8  and  4.9  show  the  pres¬ 
sure  evolution  in  the  VascoMax  300  slipper  at  1,200  m/s  and  1,500  m/s,  respectively. 
Although  the  simulation  implies  that  equivalent  plane  strain  Hugoniot  elastic  limit 
is  influenced  by  sliding  velocity,  an  implication  that  is  easily  accepted,  the  predicted 
value  of  2.8664  GPa  is  an  approximate  value  representing  this  limit. 

4-4  Validation  of  Plane  Strain  Elastic  Wave  Speed 

An  equation  to  determine  the  clastic  wave  speed  through  a  solid  material  under 
plane  strain  conditions  was  provided  in  Section  2.5.2.  Equation  4.5,  derived  from 
equations  presented  by  Zukas  [37]  and  Saada  [32],  is  a  function  of  Poisson’s  ratio,  u 
and  clastic  modulus,  E. 
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Materials  at  6.60e-ll  seconds 
0.0215  I  ■  ■  i  ■  ■  ■  ■  i  ■  ■  ■  ■  i 

0.021  - 


0.069  0.07  0.071 

X  (cm) 

1080  Steel  Rail 
YascoMax  300  Slipper 

Figure  4.6:  Location  in  Model  Where  Pressure  Data  is 

Recorded 
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Pressure  at  (x,y)  =  (0.0692,0.0202) 


time  (s) 


Figure  4.7:  Pressure  Evolution  in  VascoMax  300  Slipper  at 

1,000  m/s  Sliding  Velocity 


Pressure  at  (x,y)  =  (0.0692,0.0202) 


Figure  4.8:  Pressure  Evolution  in  VascoMax  300  Slipper  at 

1,200  m/s  Sliding  Velocity 
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Pressure  at  (x,y)  =  (0.0692,0.0202) 


Figure  4.9:  Pressure  Evolution  in  VascoMax  300  Slipper  at 

1,500  m/s  Sliding  Velocity 


ce,ps  —  \  X — n  ^  9  vi  (4-5) 

V  3p0(l-2i/)(l  +  z/) 

Solving  the  equation  for  VascoMax  300,  with  Poisson’s  ratio  equal  to  0.283  and 
elastic  modulus  of  180.7  GPa,  gives  an  elastic  wave  speed  under  plane  strain  conditions 
approximately  6,230  m/s.  This  can  be  validated  by  tracing  a  pressure  wave  through 
the  material  with  respect  to  distance  and  time.  This  is  done  by  plotting  the  pressure 
along  a  diagonal  from  the  point  of  impact  with  respect  to  distance  at  each  time  step. 
A  MATLAB  code  was  created  to  do  this.  Before  the  MATLAB  code  is  used,  the  CTH 
input  deck  needed  to  be  modified.  This  modification  and  MATLAB  code  is  discussed 
in  Appendix  D.  The  MATLAB  code  calculates  and  plots  the  change  in  pressure  at  a 
point  along  the  three  diagonals  with  respect  to  time. 

Figure  4.10  shows  the  change  in  pressure  in  the  VascoMax  300  slipper  with 
respect  to  time  on  a  45°  diagonal  at  30  /rm  from  the  point  of  impact  with  a  sliding 
velocity  of  1,000  m/s.  Figure  4.11  shows  a  similar  plot  along  the  45°  diagonal  at 
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60  /mi.  These  two  plots  can  be  used  to  determine  the  speed  of  the  elastic-plastic 
wave  generated  during  the  CTH  simulation.  The  speed  of  the  wave  is  determined 
by  dividing  the  difference  in  distance  along  the  diagonals  by  the  difference  in  time 
at  which  a  constant  value  of  pressure  is  achieved  in  the  two  figures.  In  this  case, 
the  speed  of  a  wave  with  a  constant  pressure  of  -5  GPa  was  selected.  The  change 
in  distance  along  the  diagonal  is  30  /mi,  and  the  change  in  time  is  approximately 
4.90  x  10-9  seconds.  This  gives  an  elastic-plastic  wave  speed  of  6,120  m/s. 


Pressure  at  30  |jm  Along  45  Diagonal 


Figure  4.10:  Pressure  at  30  /tin  on  a  45°  diagonal  at  1,000 
m/s  Sliding  Velocity 

The  predicted  clastic  speed  wave  through  VascoMax  300  under  plane  strain 
conditions,  given  by  Equation  4.5,  is  6,230  m/s.  This  is  the  speed  of  a  purely  elastic 
wave.  The  wave  speed  calculated  from  the  CTH  simulations  represents  an  elastic- 
plastic  wave.  Introducing  plasticity  reduces  the  speed  of  the  wave.  Therefore,  an 
elastic-plastic  wave  speed  through  the  VascoMax  300  slipper  of  6,120  m/s  suggests 
that  Equation  4.5  is  a  valid  approach  to  calculating  the  speed  of  an  elastic  wave 
through  a  solid  material. 
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Pressure  at  60  |jm  Along  45  Diagonal 


Figure  4.11:  Pressure  at  60  /mi  on  a  45°  diagonal  at  1,000 
m/s  Sliding  Velocity 

4-5  Equation  of  State  at  Low  Velocities 

Section  3.2.4  discussed  modeling  issues  when  using  an  equation  of  state  in  low 
velocity  impact  problems.  Numerical  inconsistencies  were  observed  when  solving  for 
the  mechanical  wear  rates  at  50  m/s  and  100  m/s.  These  wear  rate  values  were  two  to 
three  times  greater  than  simulations  at  higher  velocities.  The  low  velocity  simulations 
also  produced  an  inconsistent  state  of  pressure  in  the  materials.  When  the  slipper 
collides  with  the  asperity,  a  pressure  wave  is  created.  Figure  4.12  shows  a  pressure 
wave  generated  by  the  collision  of  the  VasocoMax  300  slipper  sliding  at  1,500  m/s  into 
the  surface  asperity.  This  pressure  wave  extends  into  both  the  VascoMax  300  slipper 
and  the  AISI  1080  steel  rail.  Figure  4.13  shows  the  inconsistent  state  of  pressure  at 
sliding  velocity  of  50  m/s.  Since  the  equation  of  state  is  used  to  solve  for  pressure,  it 
can  be  assumed  that  the  inconsistencies  are  a  result  of  an  improper  EOS  model. 

An  EOS  must  be  defined  for  each  material  in  CTH.  Zukas  [37]  suggests  a  mod¬ 
ification  to  the  Mie-Gruneisen  for  low  velocity  impact.  This  modification,  shown  in 
Equation  4.6),  sets  the  pressure  to  a  product  of  the  bulk  modulus,  K,  and  plastic 
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Figure  4.12:  Pressure  Wave  Generated  by  1,500  m/s  Collision 


X  (cm) 

Figure  4.13:  Inconsistent  State  of  Pressure  at  50  m/s  Sliding 
Velocity 
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strain,  ji.  It  is  reasonable  to  expect  that  a  low  velocity  impact  will  produce  less  de¬ 
formation  in  the  material.  The  proposed  equation  is  used  to  ensure  that  a  state  of 
zero  pressure  is  achieved  with  zero  compression. 

P  =  K  /i  (4.6) 

Based  on  a  review  of  the  CTH  user’s  manual  [18]  and  discussions  with  the 
author,  an  approach  was  developed  to  use  the  modified  Mie-Griineisen  EOS  by  altering 
user  input  variables  for  the  existing  CTH  Mie-Griineisen  model.  A  user  defined  Mie- 
Griineisen  EOS  in  CTH  must  define  the  material  density,  po,  sound  speed  through 
the  material,  cs,  a  linear  coefficient  in  the  Hugoniot  fit,  S,  the  Griineisen  constant, 
T,  and  the  specific  heat,  cv.  The  material  density  and  sound  speed  do  not  change 
for  the  proposed  modification.  However,  if  the  S  value  is  set  to  a  small  value,  but 
not  zero,  the  sound  speed  dependence  on  pressure  is  removed.  Also,  if  the  Griineisen 
constant  is  set  to  a  small  number,  but  not  zero,  and  the  specific  heat  is  set  to  a 
large  number,  the  material  should  be  prevented  from  changing  temperature  during 
compression  or  expansion.  This  results  in  a  model  that  should  allow  the  material  to 
respond  elastically  with  respect  to  its  bulk  modulus  and  remove  the  thermal  portion 
of  the  EOS.  This  attempt  at  implementing  the  modified  Mie-Griineisen  EOS  did  not 
remove  the  numerical  inconsistencies  at  low  velocities. 

It  was  decided  that  the  Mie-Griineisen  EOS  is  not  adequate  for  the  model  devel¬ 
oped  in  this  research.  The  Sesame  equation  of  state,  discussed  in  Section  2. 8. 2. 2  was 
selected  for  use  in  this  research  because  the  model  uses  experimental  data.  It  should 
be  noted  that  the  Sesame  EOS  interpolates  data  between  experimental  data  when 
the  sliding  velocity  falls  within  two  points.  When  the  sliding  velocity  is  outside  the 
range  of  experimental  data,  the  Sesame  EOS  extrapolates  using  experimental  data. 
The  accuracy  of  the  model  at  varying  velocity  is  dependent  on  the  amount  of  exper¬ 
imental  data.  The  lower  velocities  are  certainly  extrapolated  from  the  experimental 
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data,  and  therefore  some  error  may  be  introduced  when  running  simulations  outside 
of  the  experimental  range. 

4-6  Wolf  son  Data 

As  discussed  in  Section  1.3,  Wolfson  [36]  ran  experiments  to  study  the  wear  of 
materials  in  high  speed  track  applications.  Of  the  sixty  tests,  the  results  of  two  can 
be  used  to  compare  against  the  model  presented  in  this  thesis.  The  two  tests  both  use 
specimens  made  of  stainless  steel,  and  a  bare  steel  track  with  welded  joints.  This  is  a 
close  representation  of  the  VascoMax  300  on  AISI  1080  steel  sliding  scenario  studied 
in  this  research.  Table  4.3  shows  the  results  of  the  two  experiments.  The  experimental 
average  wear  rates  are  given  in  units  of  in/ft.  This  was  measured  on  specimens  with 
a  constant  contact  area,  An,  of  1  square  inch  (645.16  mm2).  Equation  4.7  is  used  to 
convert  the  experimental  average  wear  rates  to  units  of  mm3 /mm. 


II  ^VwolfsonA-n 


(4.7) 


Table  4.3:  Data.  From  Wolfson’s  Experiments  [36] 


Sliding 

Velocity 

(ft/s) 

Average 
Wear  Rate 
(in/ft) 

Average 
Wear  Rate 
(mm3  /mm) 

825 

2,500 

9.50  x  10"6 

7.50  x  10"6 

5.11  x  10~4 
4.03  x  10”4 

It  is  important  to  note  that  Wolfson’s  experiments  produce  three  dimensional 
wear  rates.  Therefore,  a  conversion  method  must  be  applied  to  better  represent  a 
plane  strain  environment.  To  do  this,  Archard’s  equation  [3,4],  Equation  4.8  is  used, 
which  solves  for  the  volume  of  worn  material  at  low  velocities,  W a-  hi  this  equation, 
P  is  the  applied  normal  pressure,  An  is  the  contact  area,  Jza  is  the  dimensionless 
Archard  wear  coefficient,  and  H  is  the  material  hardness.  An  Archard  wear  coefficient 
of  4.40  xl0~5  is  used  for  low  speed  wear,  and  the  VascoMax  300  material  hardness  is 
0.5  xlO3  Pa. 
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(4.8) 


WA  = 


When  solving  for  Equation  4.8  using  data  from  Wolfson’s  experiments,  which 
applied  a  constant  normal  pressure  of  300  psi,  using  a  pin  with  a  constant  contact 
area  of  1  square  inch,  a  wear  rate  value  of  1.17xl0-4  nun2  is  found.  Results  from 
Hale’s  dissertation  [20]  at  10  m/s  give  an  area  of  worn  material  equal  to  2.65  x  10~5 
nun2.  Dividing  Wolfson’s  worn  area  by  Hale’s  provides  a  constant,  Nwoifson,  which 
is  used  to  relate  Wolfson’s  experimental  data  [36]  to  the  plane  strain  model  using 
Equation  4.9. 


A^ps^wolfson  at  (4.9) 

4^  wolf  son 

where  Wps, wolf  son  is  the  plane  strain  equivalent  of  Wolfson’s  experimental  data  and 
W  is  the  Wolfson’s  experimental  wear  data  converted  from  English  to  metric  units. 
The  two  experiments  were  run  at  different  velocities  (252  m/s  and  762  m/s).  This 
provides  two  data  points  representing  Wolfson’s  experiments  converted  to  plane  strain 
wear.  These  two  data  points  are  provided  in  Table  4.4  and  will  be  plotted  along  with 
estimated  wear  rates  from  the  CTH  simulation  in  Section  4.7. 


Tabic  4.4:  Wolfson  Plane  Strain  Wear  Rates 


Sliding  Velocity 

Mechanical  Wear  Rate 

(m/s) 

(mm3  /mm) 

252 

1.16  x  1(T4 

762 

9.13  x  10~5 

Jf.,1  Mechanical  Wear  Rate  Results 

Mechanical  wear  rates  are  calculated  by  passing  output  data  from  CTH  into 
the  MATLAB  code  described  in  Appendix  C.  Figure  4.14  is  a  plot  of  the  calculated 
mechanical  wear  rates  due  to  the  collision  with  a  single  semi-circular  surface  asperity 
in  plane  strain  conditions.  The  figure  shows  four  lines  representing  the  two  failure 
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criteria,  strain  at  max  stress  and  critical  von  Mises  stress,  with  a  ”  slide”  line  interface 
condition,  and  a  “no  slide”  condition.  The  Wolfson  data  [36],  converted  to  the  plane 
strain  scenario  in  Section  4.6,  is  included  in  Figure  4.14  along  with  the  wear  rate 
results  from  Hale’s  FEA  model  [20] .  Including  the  Wolfson  data  provides  a  validation 
to  the  simulation  results  because  it  shows  the  calculated  values  are  on  the  same  order 
of  magnitude  as  an  experimental  set  of  data.  The  wear  rate  data  are  given  in  Table  4.5. 

There  are  two  failure  criteria  and  two  interface  boundary  conditions  evaluated  in 
this  thesis.  The  boundary  condition  has  an  obvious  effect  on  the  calculated  wear  rates. 
The  “slide  line”  condition  simulates  a  frictionless  surface  by  setting  the  shear  stress 
along  the  surface  to  zero,  whereas  the  “no  slide”  condition  represents  a  surface  with 
a  semi-infinite  coefficient  of  friction.  The  “no  slide”  boundary  condition  estimates  a 
higher  wear  rate  than  the  “slide  line”  condition  for  both  failure  criteria.  This  is  due  to 
the  fact  that  the  “no  slide”  condition  requires  pressure  along  the  interface  to  reach  a 
threshold  before  the  material  can  move.  Some  of  this  additional  pressure  is  captured 
within  the  data  collection  and  adds  damaged  material  to  the  calculation  resulting  in 
a  higher  wear  rate. 

Regardless  of  the  interface  boundary  condition,  the  critical  von  Mises  stress 
criteria  tends  to  result  in  a  higher  calculated  mechanical  wear  rate  than  the  strain  at 
max  stress  criteria.  The  critical  von  Mises  stress  failure  criteria  causes  the  wear  rate 
to  increase  with  sliding  velocity  up  to  1,300  m/s,  then  a  slight  decrease  is  observed. 
Wear  rates  calculated  using  the  strain  at  max  stress  failure  criteria  appear  to  level  off 
at  velocities  above  1,200  m/s.  Both  failure  criteria  follow  a  similar  curve  provided  by 
Hale  until  the  wear  rate  calculated  by  the  FEA  approach  decreases  with  increasing 
velocity  above  600  m/s.  Although  each  curve  provides  a  different  estimation  for 
mechanical  wear  rates  based  on  velocity,  they  all  are  on  the  same  order  of  magnitude. 
The  total  mechanical  wear,  presented  in  Section  4.8,  provides  a  better  overall  analysis 
of  the  models  presented  in  this  thesis. 
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■  Slide  Line  Strain  at  omax 
-IHSlide  Line  von  Mises  stress 
— *^No  slide  Strain  at  omax 
■  *  No  Slide  von  Mises  stress 
-4— Hale(FEA) 

•  Wolfson  Data 


Velocity  (m/s) 


Figure  4.14:  Plane  Strain  Mechanical  Wear  Rates 


Table  4.5:  Tabulated  Wear  Rates 


Horizontal 

Sliding 

Velocity 

(m/s) 

No  Slide 
Strain  at 
Max  Stress 
(mm3/ mm) 

Slide  Line 
Strain  at 
Max  Stress 
(mm3/ mm) 

No  Slide 
Critical  Von 
Mises  Stress 
(mm3 /mm) 

Slide  Line 
Critical  Von 
Mises  Stress 
(mm3 /mm) 

200 

1.51  x  10"4 

7.49  x  10"5 

1.26  x  10~4 

3.71  x  HT5 

300 

1.76  x  10"4 

8.34  x  10"5 

1.97  x  HT4 

6.16  x  10”5 

400 

1.93  x  10"4 

1.05  x  10"4 

2.55  x  HT4 

1.09  x  HT4 

500 

2.23  x  10"4 

1.45  x  10"4 

3.19  x  HT4 

1.70  x  HT4 

600 

2.43  x  10"4 

1.70  x  10"4 

3.49  x  10~4 

2.04  x  10~4 

700 

2.61  x  10"4 

1.86  x  HT4 

3.52  x  10'4 

2.42  x  10~4 

800 

2.66  x  10"4 

2.06  x  HT4 

3.54  x  10-4 

2.79  x  10~4 

900 

2.71  x  10"4 

2.20  x  HT4 

3.57  x  10-4 

3.03  x  lO’4 

1,000 

2.76  x  10"4 

2.32  x  HT4 

3.66  x  HT4 

3.21  x  10-4 

1,100 

2.79  x  10"4 

2.44  x  10~4 

3.72  x  10"4 

3.33  x  lO’4 

1,200 

2.86  x  10"4 

2.58  x  10~4 

3.73  x  10'4 

3.44  x  10~4 

1,300 

2.91  x  10"4 

2.62  x  10~4 

3.80  x  10”4 

3.63  x  10~4 

1,400 

2.89  x  10"4 

2.64  x  10~4 

3.77  x  10"4 

3.61  x  10~4 

1,500 

2.88  x  10"4 

2.67  x  10~4 

3.67  x  10”4 

3.45  x  10~4 
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4-8  Total  Mechanical  Wear  Results 


As  mentioned  in  Section  3.4,  the  total  mechanical  wear  of  an  HHSTT  slipper 
can  be  determined  by  plotting  the  wear  rates  as  a  function  of  distance  along  the 
track  and  integrating  with  respect  to  distance.  Figure  4.15  shows  the  estimated 
total  mechanical  wear  for  the  four  cases  presented  in  Section  4.7  as  well  as  the  total 
experimental  wear  from  the  January  2008  test  mission.  The  experimental  wear  was 
determined  by  measuring  the  thickness  of  the  slipper  at  the  end  of  the  test  compared 
to  the  design  nominal  thickness.  There  are  two  things  to  consider  when  using  this 
experimental  value.  The  first  is  that  the  third  sled  reaches  a  maximum  velocity  close 
to  1,500  m/s  and  then  decelerates  to  approximately  600  m/s  (1,342  miles  per  hour) 
at  the  end  of  the  track,  at  which  point  the  sled  and  slippers  leave  the  track,  bouncing 
along  the  ground.  It  is  not  unreasonable  to  assume  that  some  wear  occurs  during  this 
time.  The  second  thing  to  consider  is  that  the  slippers  are  not  measured  prior  to  the 
test.  The  initial  thickness  of  the  slippers  is  assumed  to  be  the  nominal  thickness  from 
the  HHSTT  design  manual  [1].  Since  the  total  volume  of  worn  material  is  determined 
as  units  of  mm3,  a  slight  deviation  from  the  nominal  thickness  can  have  an  effect  on 
the  experimental  wear  value.  Hale  [20]  gives  the  total  wear  volume  from  the  aft  right 
slipper  on  the  third  stage  sled  from  the  January  2008  test  mission  as  10,516  mm3.  The 
total  volume  of  worn  material  and  percentage  of  experimental  wear  for  each  criteria 
is  given  in  Table  4.6. 


Table  4.6:  Estimated  Total  Mechanical 

Wear 

Failure  Criteria 

Volume  of  Worn 
Material  (mm3) 

Percentage  of 
Experimental  Wear 

No  Slide  Strain  at  Max  Stress 

6,418 

61.03 

Slide  Line  Strain  at  Max  Stress 

5,186 

49.31 

No  Slide  Critical  Von  Mises  Stress 

8,504 

80.87 

Slide  Line  Critical  Von  Mises  Stress 

6,857 

65.20 

Hale  FEA  method  [20] 

4,298 

40.87 
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Figure  4.15:  Total  Mechanical  Wear 


76 


4-9  Summary  of  Results 

This  chapter  discussed  the  addition  of  a  dead  load  to  the  plane  strain  model 
to  represent  the  variable  vertical  force  of  the  slipper  as  it  slides  along  the  rail.  The 
addition  of  the  dead  load  produced  wear  rates  identical  to  simulations  without  the 
dead  load.  Therefore,  the  dead  load  was  dropped  from  the  model  and  not  included 
in  the  analysis.  It  is  possible  that  the  dead  load  is  not  adequately  represented  in  the 
submodel  developed  for  this  thesis.  It  may  produce  results  if  a  larger-scale  model  was 
developed  with  an  entire  slipper  and  dead  load  attached  to  the  top. 

The  equivalent  plane  strain  Hugoniot  clastic  limit,  derived  in  Section  2.5.2,  was 
evaluated  and  compared  to  values  obtained  from  CTH  simulations.  The  equation  used 
to  determine  <Jhel,ps ,  Equation  2.15,  assumes  it  is  independent  of  sliding  velocity, 
and  gives  a  value  of  2.8664  GPa  for  VascoMax  300.  The  CTH  simulations  suggest  that 
the  onset  of  plasticity  is  influenced  by  the  sliding  velocity  of  the  slipper.  However, 
2.8664  GPa  appears  to  be  a  good  approximation  for  VascoMax  300  cthel,ps ■  This 
chapter  also  presented  a  validation  of  the  plane  strain  elastic  wave  speed,  Equation 
2.16,  derived  in  Section  2.5.2. 

This  chapter  included  a  discussion  of  issues  encountered  during  low  velocity 
simulations.  It  is  believed  that  these  issues  are  a  result  of  the  use  of  an  improper 
equation  of  state.  A  modified  EOS  was  attempted  based  on  suggestions  by  Zukas  [37]. 
However,  this  modified  EOS  did  not  provide  a  substantial  improvement  over  the 
existing  models.  Hydrocodes  are  typically  used  to  model  high  energy  problems,  such 
as  explosives  or  high  velocity  impact  problems.  Therefore,  the  low  velocity  issues 
were  not  a  significant  surprise.  The  simulations  ran  with  no  errors  over  a  velocity 
range  from  200  m/s  to  1,500  m/s.  Simulations  with  a  sliding  velocity  below  200  m/s 
resulted  in  numerical  inconsistencies  described  in  Section  4.5. 

This  chapter  also  presented  the  results  of  the  CTH  simulation  for  four  scenarios; 
“No  slide”  strain  at  max  stress,  “slide  line”  strain  at  max  stress,  “no  slide”  critical 
von  Mises  stress,  and  “slide  line”  critical  von  Mises  stress.  The  mechanical  wear  rate 
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results  are  presented  in  Figure  4.14  and  the  total  mechanical  wear  can  be  found  in 
Figure  4.15.  The  wear  rate  data  are  given  in  Table  4.5. 
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V.  Summary  and  Conclusions 


The  chapter  is  a  summation  of  the  material  presented  in  this  thesis.  The  lit¬ 
erature  search  and  theoretical  background  will  be  discussed  first,  followed  by  a  brief 
description  of  the  hydrocode  model  developed  and  the  results  obtained.  Finally,  con¬ 
clusions  will  be  presented  on  the  results  of  the  thesis,  and  suggestions  for  future  work 
will  be  presented. 

5. 1  Summary 

Research  into  the  onset  of  wear  of  sliding  bodies  has  produced  low-velocity 
models  capable  of  estimating  worn  material.  One  model  in  particular,  the  Archard 
Wear  model  [2-4],  has  been  used  to  establish  relationships  between  wear  in  plane 
strain  to  three-dimensional  wear.  Previous  work  by  Hale  [20]  and  Meador  [28]  has 
made  use  of  a  plane  strain  scenario  to  model  the  slipper-rail  sliding  event. 

Based  on  the  previous  research,  a  hydrocode  model  was  developed  using  CTH  to 
estimate  plane  strain  mechanical  wear  rates.  The  model  allows  a  VascoMax  300  slipper 
to  collide  with  a  6  /jm  radius  surface  asperity  made  of  AISI  1080  steel.  Damage  was 
recorded  per  sliding  distance  to  give  wear  rates.  Two  failure  criteria  were  evaluated 
(Section  2.9):  critical  von  Mises  stress  and  strain  at  max  stress.  These  failure  criteria 
were  established  by  the  Johnson-Cook  viscoplasticity  model  presented  in  Section  2.4. 
The  model  also  has  two  distinctly  different  interface  boundary  conditions  between 
the  slipper  and  rail.  One  boundary  condition,  “slide  line”,  simulates  a  frictionless 
surface  by  setting  the  shear  stress  along  the  surface  to  zero.  The  second  boundary 
condition,  “no  slide”,  simulates  a  surface  with  a  semi-infinite  coefficient  of  friction 
by  establishing  a  pressure  threshold  that  must  be  exceeded  for  the  material  to  move. 
These  two  boundary  conditions  represent  two  extremes  along  the  surface  in  terms  of 
friction. 

Since  the  model  developed  in  this  thesis  simulates  a  collision  between  two  metals 
at  high  velocities,  and  establishes  failure  criteria  based  on  the  material  response  due 
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to  the  collision,  fundamentals  of  wave  mechanics  were  researched.  Previous  research 
of  high  velocity  impact  has  considered  uniaxial  strain  conditions.  As  such,  equations 
exist  that  are  used  to  estimate  the  onset  of  plasticity  in  the  material,  called  the  Hugo- 
niot  elastic  limit,  and  to  calculate  the  speed  of  pressure  waves  propagating  through 
a  solid  material  under  uniaxial  strain  conditions.  However,  the  scenario  of  interest 
for  this  research  is  plane  strain.  Therefore,  equations  were  derived  in  Chapter  II  to 
evaluate  the  plane  strain  elastic  wave  speed  through  a  solid,  Equation(2.16),  and  the 
equivalent  plane  strain  Hugoniot  elastic  limit,  Equation  2.15. 

The  CTH  model  was  run  at  velocities  ranging  from  200  m/s  to  1,500  m/s.  The 
results  of  the  simulations  were  presented  in  Chapter  IV.  The  estimated  mechanical 
wear  rates  were  used  to  determine  the  total  mechanical  wear  of  the  aft  right  slipper 
from  the  third  sled  of  the  January  2008  test  mission. 

5.2  Conclusions 

The  plane  strain  derivations  presented  in  Chapter  II  were  validated  in  Chap¬ 
ter  III  using  results  obtained  from  CTH  simulations.  The  derived  equations  imply 
they  are  independent  of  sliding  velocity.  However,  the  CTH  simulations  suggest  that 
onset  of  plasticity  and  elastic  wave  speed  through  the  slipper  are  effected  by  the  slid¬ 
ing  velocity.  This  effect  appears  to  be  minimal,  and  the  derived  equations  present  an 
adequate  estimation  of  these  values. 

The  addition  of  a  dead  load  to  represent  the  vertical  force  of  the  slipper  into 
the  rail  was  presented  in  Section  4.1.  The  dead  load  was  added  to  provide  a  more 
accurate  representation  of  the  HHSTT  environment.  CTH  does  not  allow  a  force 
input.  Therefore,  Equation  4.1  was  used  to  calculate  the  effective  mass  representing 
the  vertical  force  from  the  DADS  data.  The  simulations  with  the  dead  load  produced 
mechanical  wear  rates  that  were  identical  to  simulations  without.  One  possible  rea¬ 
son  for  these  results  is  the  fact  that  a  submodel  of  the  slipper  was  developed,  only 
representing  a  small  part  of  the  slipper  in  the  simulation.  It  may  be  necessary  to 
develop  a  full  model  of  the  slipper  with  an  attached  dead  load  to  get  an  accurate 
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respresentation  of  the  vertical  force.  Such  a  model  was  not  developed  for  this  thesis, 
because  the  simulation  time  would  drastically  increase. 

The  mechanical  wear  rates  obtained  from  the  CTH  simulation  appear  to  be 
an  accurate  estimation  of  the  HHSTT  slipper-rail  sliding  event.  Data  from  Wolfson 
[36]  was  converted  to  the  plane  strain  scenario  in  Section  4.6  and  plotted  with  the 
estimated  mechanical  wear  rates  in  Figure  4.14.  The  two  data  points  are  of  the  same 
order  of  magnitude  as  the  estimated  values  from  the  CTH  simulation.  The  total 
mechanical  wear  was  calculated  and  plotted  in  Figure  4.15.  The  total  experimental 
wear  from  the  January  2008  test  mission  was  determined  by  Hale  [20]  to  be  10,516 
nun3.  The  results  of  the  total  mechanical  wear  predict  between  49.31%  and  80.87% 
of  the  experimental  wear.  However,  due  to  the  uncertainty  of  the  true  experimental 
wear,  discussed  in  Section  4.8,  the  results  of  the  simulation  are  acceptable.  The  results 
of  the  simulation  and  total  wear  calculation  suggest  that  the  model  developed  is  an 
adequate  method  to  model  mechanical  wear. 

5.3  Future  Work  Suggestions 

The  model  presented  in  this  thesis  has  been  developed  and  modified  based  on 
previous  research.  There  are  however,  simplifications  made  to  allow  the  plane  strain 
model  to  represent  a  three-dimensional  scenario.  One  significant  simplification  to 
this  research  is  the  absence  of  a  thermal  model.  Previous  work  by  Meador  [28]  has 
shown  that  melt  wear  plays  an  important  role  in  the  slipper  -  rail  sliding  event.  Mrs. 
Gracie  Paek,  as  part  of  her  PhD  dissertation  is  developing  a  thermodynamic  model 
to  represent  melt  wear  of  HHSTT  slippers. 

The  Johnson-Cook  fracture  model,  discussed  in  Section  2.9.3,  can  be  used  to 
model  material  damage.  This  failure  criteria  was  evaulated  for  use  in  the  current 
model.  However,  the  simulations  estimated  zero  wear  when  this  criteria  was  used. 
The  fracture  model  requires  five  coefficients  to  be  defined  for  each  material.  These 
fracture  coefficients  are  not  known  for  VascoMax  300  or  AISI  1080  steel.  AISI  4340 
steel  and  iron  were  used  to  represent  the  VascoMax  300  slipper  and  AISI  1080  steel 
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rail,  respectively.  With  properly  defined  fracture  coefficients,  obtained  through  ex¬ 
perimentation,  the  Johnson-Cook  fracture  model  may  yield  worn  material. 
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Appendix  A.  Plane  Strain  Derivations  of  Hug oniot  Elastic  Limit  and 

Elastic  Wave  Speed 

An  approach  for  considering  the  uniaxial  Hugoniot  elastic  limit  is  outlined  in  Meyers 
[29]  and  Zukas  [37].  A  similar  process  is  applied,  but  in  this  process  a  Hugoniot 
clastic  limit  is  determined  considering  plane  strain  conditions.  The  speed  of  the 
clastic-plastic  wave  generated  clue  to  plane  strain  collision  is  also  determined.  The 
following  equations  show  the  derivation  of  Equations  (2.15)  and  (2.16). 

A.l  Equivalent  Hugoniot  Elastic  Limit  for  Plane  Strain 

Since  plane  strain  is  considered,  one  can  say 


£1  —  e\  +  4 

(A.l) 

£2  =  4  +  4 

(A.2) 

=4+4=0 

(A.3) 

£3  =  ~4 

(A.4) 

where  e\  =  elastic  strain  and  ef  =  plastic  strain.  The  next  step  is  to  consider  the 
plastic  portion  to  be  incompressible,  thus 


£i  +  4  +  4  =  0  (A. 5) 

4  +  4  =  4  (A. 6) 

If  Equations  (A.l)  through  (A. 6)  are  combined,  the  result  is  an  equation  for  the 
summation  of  principal  strains  £\  and  £2- 
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£l  +  £2  —  £3  —  £1  —  £1  +  £2  —  £2 

(A.7) 

(£1  +  £2)  =  £1  +  £2  +  £3 

(A.8) 

Saada  [32]  provides  constitutive  equations  for  plane  strain. 

£i  =  E  [(!  v)ax  ua2] 

(A.9) 

4=  E  K1  wj 

(A. 10) 

^3  =  0 

(A. 11) 

The  Tresca  yield  theory  is  used  to  get  a2  =  f(cri,  Y0)  [32], 

Y(t  =  a1-a2 

(A. 12) 

a2=(J\-  Y0 

(A. 13) 

where  Y0  is  the  yield  stress  for  a  uniaxial  elastic  -  perfectly  plastic  material.  The 
constiutive  equations  along  with  Tresca  yield  theory  are  used  to  reduce  Equation  A. 8 
to  a  form  for  oq. 


(ei  +  £2  )E  Y0 

(1  +  u) (1  -  2u)2  +  ~2 


(A. 14) 


The  next  step  is  to  consider  loading  of  an  clastic  -  perfectly  plastic  material, 
starting  with  the  pressure. 
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(A. 15) 
(A. 16) 


p  _  01  +  02  +  03 

“  3 

p  _  01  +  02  +  Kgl  +  02) 

3 

where  z/  is  the  Poisson’s  ratio,  and  03  =  1/(01  +  02)  in  a  plane  strain  scenario.  If  one 
solves  again  for  ay, 


=  3  P  Y0 
2  (1  +  v)  +  2 


(A. 17) 


0i  becomes  the  stress  of  importance  for  the  analysis  as  it  is  in  a  uniaxial  strain 
situation.  When  pressure,  P,  equals  zero 


(A.i8) 


Figure  A.l  shows  the  loading  of  an  elastic-perfectly  plastic  material.  For  the 
case  of  zero  pressure 


(A. 19) 


For  our  purposes,  we  assume  f3  =  |/i,  as  related  to  a  uniaxial  stress  situation  [32], 
where  //.  =  yqy-yp  The  summation  of  principal  strain  term  £\  +  £2,  is  evaluated  in 
terms  of  a  and  f3  to  get 


(£1  +  £2) 


3  Y0 

4  E 


(1  +  *) 


(A. 20) 


Equations  (A. 14)  and  (A. 20)  are  combined  to  solve  for  the  equivalent  Hugoniot 
elastic  limit  for  the  case  of  plane  strain. 
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Figure  A.l:  Loading  of  an  Elastic-Perfectly  Plastic  Material 


&HEL,PS  ~  Yo 


8  -  16z/ 


(A. 21) 


A. 2  Plane  Strain  Elastic  Wave  Speed 

Combining  equations  (A. 13)  and  (A. 20),  the  result  for  (ex  +  e2)  is 


(Cl  +  £2) 


3(l  +  i/) 
4  E 


(Vl  +  02) 


(A. 22) 


and 


v 

cr2  =  - - a 

1  —  v 


(A. 23) 


Combining  equations  (A. 22)  and(A.23),  and  solving  for  a±,  the  following  is 


obtained 


(A. 24) 


cr  i 


4(1  ~  v) 

3(1  +  v)(l  —  2u) 


(el  +  £2)E 


The  speed  of  sound  through  any  medium  can  be  represented  as  the  square  root 
of  pressure  divided  by  density.  Through  a  solid  media,  the  pressure  term  is  replaced 
by  the  clastic  modulus.  For  the  case  of  plane  strain,  the  clastic  modulus  carries  an 
added  term  representing  the  summation  of  strain,  (eq  +  £2).  The  elastic  wave  speed 
for  the  case  of  plane  strain  is  given  by  Equation  A. 25 


CE,PS 


4  (1  ~  u) 

3  Po(l  -  2u)(l  +  u) 


E 


(A. 25) 
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Appendix  B.  CTH  Input  Deck 


B.l  Discussion  of  CTH  Input 

The  sample  CTH  input  deck  provided  in  this  appendix  is  used  to  estimate  the 
wear  rate  of  VascoMax  300  sliding  at  1,000  m/s  into  a  6  pm  radius  surface  asperity 
made  of  AISI  1080  steel.  This  input  deck  can  be  altered  to  represent  an  event  with  a 
different  sliding  velocity.  To  do  this,  the  total  simulation  time,  “tstop” ,  and  intermedi¬ 
ate  step  time  need  to  be  updated.  The  “tstop”  variable  is  found  on  line  36  of  the  input 
deck.  There  are  nine  variables  in  the  input  deck  that  use  the  intermediate  step  time, 
and  all  nine  must  be  updated  for  each  run.  The  variables  are:  “dt” ,  “dtfrequency” , 
“PlotTime”,  “SaveTime”,  and  “HisTime”.  The  “dt”  variable  is  listed  four  times  on 
lines  244,  247,  318,  and  322.  The  “dtfrequency”  variable  is  listed  twice  on  lines  250 
and  253.  The  “PlotTime”  variable  is  found  on  line  330,  the  “SaveTime”  variable  is 
found  on  line  331,  and  the  “HisTime”  variable  is  found  on  line  406. 

The  initial  velocity  also  needs  to  be  updated  when  changing  the  sliding  velocity 
of  the  VascoMax  300  slipper.  The  initial  velocity  vector  is  found  in  the  diatom 
input  set,  starting  on  line  124.  Since  the  VascoMax  300  slipper  is  the  only  moving 
body  for  this  analysis,  it  is  the  only  velocity  that  needs  to  be  changed.  The  initial 
velocity  input  for  the  VascoMax  300  slipper  is  found  on  line  152  of  the  input  deck. 
The  initial  velocities  of  the  other  materials  should  remain  at  zero  for  all  simulations. 
It  is  important  to  note  that  the  base  unit  of  length  in  CTH  is  cm.  This  means 
that  all  velocities  and  distances  must  be  input  having  dimensions  of  cm/s  and  cm, 
respectively.  This  is  seen  on  line  152  of  the  included  input  deck,  the  initial  velocity 
is  given  as  100,000  cm/s  which  gives  1,000  m/s. 

In  CTH,  the  units  of  pressure  and  stress  are  expressed  as  dynes/cm2  and  tem¬ 
perature  in  electron  volts,  eV .  This  is  accounted  for  in  the  post-processing  code 
provided  in  Appendix  C.  However,  the  implementation  of  the  Johnson-Cook  vis¬ 
coplastic  model  requires  the  conversion  of  a  couple  material  constants.  Pressure  and 
stress  is  converted  to  GPa  from  dynes/cm2  using  Equation  (B.l),  and  temperature 


is  converted  from  K  to  eV  using  Equation  (B.2)  [23].  Table  B.l  shows  the  material 
constants  in  units  compatible  with  CTH. 


F dynes/ cm?  =  PgPcl  X  1010  (B.l) 

TeV  =  Tk/  11604.505  (B.2) 


Table  B.l:  Johnson-Cook  Coefficients  for  VascoMax 
300  and  AISI  1080  Steel  in  CTH  Units  [15,20,28] 


Coefficient 

VascoMax  300 

AISI  1080  Steel 

A  {dynes /cm2) 
B  {dynes / cm2) 

C  (Unitless) 
m  (Unitless) 
n  (Unitless) 

2.1  x  101U 
0.124  x  1010 
0.03 

0.8 

0.3737 

0.7  x  101U 

3.6  x  1010 

0.17 

0.25 

0.6 

B.2  Example  CTH  Input  Deck 


*eor*  cthin 

* 

*  cthin  input  with  Spymaster  graphics  for  slipper  wear  simulation 
6  * 

*  filename:  slipperwear . in 

* 

*  1.  File  modified  by  Steve  Meador  (MS-10M) 

*  2.  File  converted  to  CTH  v8 . 1  by  Maj  Chad  Hale,  PhD-09S,  Aug  ... 

2008 

11  *  3.  new  format  based  on  CTH  Course  (4-7  Aug  08)  in  Albuquerque,  ... 
NM 

*  4.  modifies  Cameron’s  393  m/s,  No  Coating,  Asperity,  T=297  input... 

file 


* 

* 

*  |  - >  | 

16  *  I  I  I 

*  I  v  / 

*  - 

* 

* 


V 


21 

26 

31 

36 

41 

46 

51 

56 

61 

66 


*  vx=varies ,  vy=-l  m/s  V 300  Steel  Slider,  1080  Steel  Rail,  No 
Atm  . 

*  No  Slide  line.  mix=l  frac=l  Rounded  corner. 

*  Added  mass  on  top  to  simulate  sled  mass 

*  title  record  set 

Horizontal  Velocity  =  1000  m/s,  Vertical  Velocity  =  -0.50  m/s 


*  control  input  set 


control 

mmp3  *  enable  multiple  material  temperatures  and 

pressures  in  each  cell 

tstop  =  6.60e-9  *  stopping  criteria  for  time  level  -  this  . 

is  total  simulation  time 

nscycle  =  100000  *  maximum  number  of  cycles  to  be  run 

*  rdumpf  =  3600.  *  time  for  back-ups  of  restart  file  ... 

updates 


tbad  =  le30 
*  dtcourant  =  0.6 
ygravity  =  -980 
endcontrol 


*  maximum  number  of  thermodynamics  warnings 
*  Courant  condition  multiplier 

*  Acceleration  due  to  gravity  =  -9.80  m/s~2 


*  mesh  input  set 

*  geom=2DR ( rectangular  x,y) 

*  geom=2DC ( cylindrical  x=radius ,  y=axis) 

*  geom=3DR (rectangular  x,y,z) 

*  type=e  (Eulerian)  now  the  default  (CTHv8.1) 

*  x#= coordinate  range  for  plot 

*  y#= coordinate  range  for  plot 

*  dxf=width  of  first  cell  in  the  region 

*  dxl=width  of  last  cell  in  the  region 

*  n=number  of  cells  added  in  this  region 

*  w=total  width  of  this  region  in  centimeters 

*  r=ratio  of  adjacent  cell  widths 


mesh 

block  1  geom=2dr  *  coordinates  for  2D  rectangular 

Eulerian  mesh 
xO  =  0.0000 

xl  w  =  850e-4  dxf  =  1.0e-4  dxl  =  1.0e-4 
endx 

yO  =  0.0000 

yl  w  =  850e-4  dyf  =  1.0e-4  dyl  =  1.0e-4 


90 


endy 

endblock 

endmesh 
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*  EOS  input  set 
eos 


76  materiall  ses  grepxyl 
Cameron ) 

material2  ses  iron 
*  MAT3  MGRUN  =  user  R0=8.13  CS 
modified  Mei -Gruneisen  for 
material3  ses  steel_v300 
material4  mgr  platinum 
81  endeos 


*  epoxy  rail  coating  (Cinnamon/... 

*  1080  steel  rail 

3 . 63  e5  SI  =  1 e -3  G0  =  le-3  CV  =  lel5  *  . 
VM300  slipper 

*  VascoMax  300  slipper 

*  platinum  for  simulated  sled  mass 


*  elastic -plastic  input  set 
86  epdata 

vpsave  *  cell  yield  stress  and  plastic  strain  rate  data  . 

is  saved 

lstrain  *  compute  and  save  Lagrangian  strain  tensor  ... 

components 

mix  =3  *  volume  averaged  yield  strength  normalized  by  sum 

of  volume  fractions 


*Epoxy  Glider  Coating 


91  matep  =  1 

poisson  0.46 
yield  1 . 0  e8 

matep  =2  *  1080 

96  JO  USER 

A JO  0 . 7el0 
BJ0  3 . 6el0 
CJ0  0.17 
M JO  0.25 
101  N  JO  0.6 

T  JO  0.14391 
poisson  0.27 


Steel  rail 

*  A 

*  B 

*  C 

*  m 

*  n 

*  Melting  temperature 


106 


111 


matep  =3  *  VascoMax  300  slipper 


JO  USER 

A  JO  =  2  .  1  e  10  * 
BJ0  =  0 . 124el0  * 
CJ0  =  0.03  * 
M JO  =  0.8  * 
N  JO  =  0.3737  * 
T  JO  =  0.145202  * 
poisson  0 . 283 


A 

B 

C 

m 

n 

Melting  temperature 
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116 

121 

126 

131 

136 

141 

146 

151 

156 

161 

166 


*  platinum  simulated  sled  mass 


matep  =  4 

poisson  0.2 
yield  lOelO 

*  SLI  2  3 

endepdata 

*  diatom  input  set 

diatom 
block  1 

package  ’1080  steel  rail’ 
material  2 
numsub  100 

temperature  =  2.55935e-2  *  eV  =  74.93F  =  297  K 
velocity  0.0,  0.0 

insert  box 
pi  0  0 

p2  850e -4  200e-4 
endinsert 
delete  circle 

center  700e-4  200e-4 
radius  6e-4 
enddelete 
insert  circle 

center  700e-4  200e-4 
radius  6e-4 
endinsert 
endpackage 

package  ’slipper’ 
material  3 
numsub  100 

temperature  =  0.0184558 
velocity  =  1000e2,  -0.50e2 

insert  box 

pi  0.0  200  e -4 
p2  694e -4  325e-4 
endinsert 
delete  box 

pi  692e -4  200e -4 
p2  694e -4  202e-4 
enddelete 
delete  circle 

center  692e-4  202e-4 
radius  2e-4 
enddelete 
insert  circle 

center  692e-4  202e-4 
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171 

176 

181 

186 

191 

196 

201 

206 

211 

216 


radius  2e-4 
endinsert 
endpackage 

endblock 

enddiatom 

*  tracer  input  set 
tracer 


add 

0 . 06755  , 

0 .01905 

to 

0.07115  , 

0 .01905 

n  =  37 

add 

0 . 06755  , 

0 .01915 

to 

0.07115  , 

0 .01915 

n  =  37 

add 

0 . 06755  , 

0 .01925 

to 

0.07115  , 

0 .01925 

n  =  37 

add 

0 . 06755  , 

0 .01935 

to 

0.07115  , 

0 .01935 

n  =  37 

add 

0 . 06755  , 

0 .01945 

to 

0.07115  , 

0 .01945 

n  =  37 

add 

0 . 06755  , 

0 .01955 

to 

0 . 07115  , 

0 .01955 

n  =  37 

add 

0 . 06755  , 

0 .01965 

to 

0 . 07115  , 

0 .01965 

n  =  37 

add 

0 . 06755  , 

0 .01975 

to 

0 . 07115  , 

0 .01975 

n  =  37 

add 

0 . 06755  , 

0 .01985 

to 

0 . 07115  , 

0 .01985 

n  =  37 

add 

0 . 06755  , 

0 .01995 

to 

0 . 07115  , 

0 .01995 

n  =  37 

add 

0 . 06755  , 

0 . 02005 

to 

0.07115  , 

0 . 02005 

n  =  37 

add 

0 . 06755  , 

0 . 02015 

to 

0.07115  , 

0 . 02015 

n  =  37 

add 

0 . 06755  , 

0 . 02025 

to 

0.07115  , 

0 . 02025 

n  =  37 

add 

0 . 06755  , 

0 . 02035 

to 

0.07115  , 

0 . 02035 

n  =  37 

add 

0 . 06755  , 

0 . 02045 

to 

0.07115  , 

0 . 02045 

n  =  37 

add 

0 . 06755  , 

0 . 02055 

to 

0 . 07115  , 

0 . 02055 

n  =  37 

add 

0 . 06755  , 

0 . 02065 

to 

0 . 07115  , 

0 . 02065 

n  =  37 

add 

0 . 06755  , 

0 . 02075 

to 

0 . 07115  , 

0 . 02075 

n  =  37 

add 

0 . 06755  , 

0 . 02085 

to 

0 . 07115  , 

0 . 02085 

n  =  37 

add 

0 . 06755  , 

0 . 02095 

to 

0 . 07115  , 

0 . 02095 

n  =  37 

add 

0 . 06755  , 

0 . 02105 

to 

0 . 07115  , 

0 . 02105 

n  =  37 

add 

0 . 06755  , 

0.02115 

to 

0.07115  , 

0.02115 

n  =  37 

add 

0 . 06755  , 

0 . 02125 

to 

0.07115  , 

0 . 02125 

n  =  37 

add 

0 . 06755  , 

0 . 02135 

to 

0.07115  , 

0 . 02135 

n  =  37 

add 

0 . 06755  , 

0 . 02145 

to 

0.07115  , 

0 . 02145 

n  =  37 

add 

0 . 06755  , 

0 . 02155 

to 

0.07115  , 

0 . 02155 

P 

II 

OJ 

->1 

add 

0 . 06755  , 

0 . 02165 

to 

0 . 07115  , 

0 . 02165 

n  =  37 

add 

0 . 06755  , 

0 . 02175 

to 

0 . 07115  , 

0 . 02175 

P 

II 

OJ 

add 

0 . 06755  , 

0 . 02185 

to 

0 . 07115  , 

0 . 02185 

n  =  37 

add 

0 . 06755  , 

0 . 02195 

to 

0 . 07115  , 

0 . 02195 

n  =  37 

add 

0 . 06755  , 

0 . 02205 

to 

0 . 07115  , 

0 . 02205 

n  =  37 

add 

0 . 06755  , 

0 . 02215 

to 

0.07115  , 

0 . 02215 

n  =  37 

add 

0 . 06755  , 

0 . 02225 

to 

0.07115  , 

0 . 02225 

n  =  37 

add 

0 . 06755  , 

0 . 02235 

to 

0.07115  , 

0 . 02235 

n  =  37 

add 

0 . 06755  , 

0 . 02245 

to 

0.07115  , 

0 . 02245 

n  =  37 

endtracer 

*  convection  control  input  set 
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221 

226 

231 

236 

241 

246 

251 

256 

261 

266 


Convct  *  enable  convection  of  internal  energy 

convection  =1  *  use  slope  of  internal  energy  and  mass 

density,  discard  KE  residual 
interface  =  smyra  *  scheme  for  interface  tracker 
endconvct 


*  fracture  input  set 

Fracts  *  enable  fracture  data  (dynes/cm~2) 

pressure 

pfracl  =  -1.0e8  *  fracture  stress  or  pressure  for  nth  . 

material 

pf rac2  =  -2 . OelO 
pfrac3  =  -7.45el0 
pf rac4  =  -1 . 2el0 

pfmix  =  -1.20el0  *  fracture  stress  or  pressure  in  a  cell 

with  no  void  present 

pfvoid  =  -1.20el0  *  fracture  stress  or  pressure  in  a  cell 

with  a  void  present 
endf  r act  s 

*  edits  input  set 

edit 
exact 
shortta 

t ime  =  0.0 
ends 
longt 

t ime  =  0 . 0  eO 
endl 
plott 

t ime  0 . 0  e -6 
endp 
histt 

t  ime  0 . 0  e -6 
htracer  all 
endhi stt 
ende 

*  boundary  condition  input  set 

*  0=symmetry 

*  l=sound  speed  based  absorbing 

*  2=extrapolated  pressure  with  no  mass  allowed  to  enter 

*  3=extrapolated  pressure  but  mass  is  allowed  to  enter 


*  short  edits  based  on  time 
dt  =  6 . 60  e - 1 1 

*  long  edits  based  on  time 
,  dt  =  6 . 60e - 1 1 

*  plot  dumps  based  on  time 
dtfrequency  6.60e-ll 

*  tracer  history  based  on  time 
dtfrequency  6.60e-ll 
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271 

276 

281 

286 

291 

296 

301 

306 

311 

316 


boundary  *  enable  boundary  condition  data 

bbydro  *  enable  hydrodynamic  boundary  . 

conditions 
block  1 

bxbot  =  1  ,  bxtop  =  2 
bybot  =  1  ,  bytop  =  2 
endb 
endh 
endb 

*  CSH :  cleaned  up  to  here... 


*heatconduction  *  enable  heat  conduction 

*  MAT1  TABLE  =3  *  conductivity  tables  defined  in 

DEFTABLE  list  below 

*  MAT2  TABLE  =  1 

*  MAT3  TABLE  =  2 

*  endh 


*  DEFTABLE=1  *  1080  STEEL 


*T(eV)  k ( erg/ s / eV/ cm ) 


*  1 . 4684e -3 

*  1 . 0377e -2 

*  1 . 9090e-2 

*  2 . 7900e-2 

*  3 . 67 1 1 e -2 

*  4.5521e-2 

*  5 . 4332  e -2 

*  6 . 3142e-2 

*  7 . 1953e-2 

*  8 . 9574e-2 

*  l.lllle-1 


4 . 7700  e 10 
4. 8100el0 
4 . 5200  e 10 
4. 1300  e 10 
3 . 8100el0 
3 . 5100el0 
3 . 2700  e 10 
3  .  OlOOelO 
2 . 4400  e  10 
2 . 6800  e 10 
3 .OlOOelO 


*  endd 


*  DEFTABLE=2  *  VascoMax  300  Steel 

*T(eV)  k ( erg/ s / eV/ cm ) 

*  3.6711e-3  2 . 4715el0 

*  1 . 4684e -2  2.7424el0 

*  2 . 9369  e -2  2.9794el0 

*  3 . 9158e-2  3.0132el0 

*  endd 

*  DEFTABLE=3  *  Epoxy 

*T(eV)  k ( erg/ s / eV/ cm ) 

*  3.6711e-3  6 . 5  e8 

*  1 . 4684e -2  6.5e8 

*  2 . 9369  e -2  6.5e8 

*  3 . 9158e-2  6.5e8 

*  endd 
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321 

326 

331 

336 

341 

346 

351 

356 

361 

366 


*mindt 

*  t ime  =  0.0  dt 

*  endm 


*  minimum  allowable  time  step  in  mesh 
6 . 60e - 1 1 


maxdt 
t  ime 
endm 


*  maximum  allowable  time  step  in  mesh 
0.0  dt  =  6 . 60e  - 1 1 


*  CSH :  Attempt  to  get  data  for  Spymaster 

spy 

PlotTime  (0 . 0 ,  6.60e-ll); 

SaveTime (0 . 0 ,  6.60e-ll); 

Save  ("  VOID  ,  V0LM  ,M  ,P  ,  XXDEV  ,  YYDEV  ,  XYDEV  ,  VX  ,  VY  ,  T  ,  TK  ,  PM  ,  TM  ,  YLD  ,  Q3  ,  J2P  .  .  . 
")  ; 

def ine  main  () 

{ 

"/,  pprintf ("  PLOT:  Cycle=%d , 

X  XLimits  (400e~4 ,725e~4) ; 

“/  YLimits  (175e~4 ,300e~4)  ; 

y  Image  (" Materials  ")  ; 
c/  Window  (0 ,0 ,0 .75 , 1)  ; 

y  Label ( sprint f ("Materials 

y  Plot 2DMa t  s  (0 . 3) ; 

‘/,  ULab  el  (  "  Test  :  (cm)"); 

y  Draw2DMesh () ; 

y  MatColors  (RED , GREEN ,  YELLOW , NO _COLOR) ; 

°/c  MatNames  ("  Epoxy  Coating  1080  Steel  Rail  "  ,  "  VascoMax  300 
Slipper  ",  "") ; 

y  DrawMatLeg end  ("",0.71,0.2,0.99,0.9); 

°/c  Endlmage; 

XLimits (650e-4,750e-4) ; 

YLimits (150e-4,250e-4) ; 

Image ("VonMisesStress  "  )  ; 

Window (0,0,0.75,1) ; 

ColorMapRange (0 ,4000)  ; 

ColorMap Clipping (OFF , OFF)  ; 

Label  ( sprintf  ("  von  Mises  Stress  at  °/,6.2e  seconds",  TIME)); 

Plot2D ( " J2P " ) ; 

Draw2DMatContour ; 

DrawColorMap ( " vonMises  Stress  (MPa)",  0 . 7 , 0 . 4 , 0 . 9  ,  0 . 9 )  ; 

Endlmage ; 

y  XLimits  (650e~4 ,750e~4) ; 

y  YLimits  (150e~4 , 250e~4) ; 

y  Imag e ( " PI  as t i cStrainRat e  " ) ; 

y  Window  (0 ,0 ,0 .75 , 1)  ; 

y  Col orMapRang  e (1 e6 , 1 el5 , LOG_MAP ) ; 


Time=y,e\n"  ,  CYCLE ,  TIME)  ; 


at  °/6.2e  seconds",  TIME)); 

°/,  toggle  on/off  mesh 
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7  Co  l  orMapCl  ipp  ing  (OFF ,  OFF)  ; 

7  Lab  el  ( sprint  f  ( "  PI  as  t  i  c  Strain  Rate  at  7,6. 2e  seconds", 

%  Plot  2D ( " PSR  ") ; 

%  Draw2DMatCont our ; 

371  7  DrawColorMap (" Plastic  Strain  Rate  (1/sec)",  0.7,0. 4,0. 

7  Endlmage; 

XLimits (685e-4,715e-4) ; 

YLimits (190e-4,215e-4) ; 

376  Image (" Materials.small  ")  ; 

Window (0,0,0.75,1) ; 

Label ( sprintf (" Mater ials  at  76. 2e  seconds",  TIME)); 
Plot2DMats  (0.3)  ; 

Label(  "Test  Label:  Distance  (cm)"  ); 

381  7  Draw2DMesh () ;  7  toggle  on/off  mesh 

MatColors (N0_C0L0R , GREEN , YELLOW , N0_C0L0R) ; 

MatNames ( " "  ,  "  1080  Steel  Rail " , " VascoMax  300  Slipper","" 
DrawMatLegend (" "  ,0.71,0.2,0.99,0.9)  ; 

Endlmage ; 

386 

XLimits (650e-4,750e-4) ; 

YLimits (150e-4,250e-4) ; 

Image ("Pressure")  ; 

Window (0,0,0.75,1)  ; 

391  ColorMap Range (le6 ,2ell , L0G_MAP ) ; 

ColorMap Clipping (OFF , OFF)  ; 

Label ( sprintf (" Pressure  at  76. 2e  seconds",  TIME)); 
Plot2D (" P")  ; 

Draw2DMatContour ; 

396  DrawColorMap (" Pressure  ( dyne / cm ~ 2 ) "  ,  0 . 7 , 0 . 4 , 0 . 9  ,  0 . 9)  ; 

Endlmage ; 

> 

401  SaveHis  ("POSITION  ,  YLD  ,  Q3  ,  PSR  ,  VOLM  +  3  ,  P  ,  XXDEV  ,  YYDEV  ,  XYDEV  ,  J2P 
7  SaveHis ("POSITION ,  YLD ,  VOLM+3 , P , XXDEV ,  YYDEV , XYDEV ") ; 

7  SaveHis (" POSITION , Q3 , PSR ,  VOLM+3") ; 

7  SaveHis ("POSITION ,  VOLM+3 , DMG3 ") ; 

SaveTracer ( ALL) ; 

406  HisTime(0,6.60e-ll); 

define  spyhi s .main ( ) 

{ 

Hi sLoad ( 1 , " hscth " ) ; 

411  LabelC'EFP  Velocity  (Tracer  1)"); 

TPlot (" VY . 1" , 1 , AUTOSCALE) ; 

> 

endspy 


TIME) ) ; 

9,0.9) ; 


)  ; 


97 


Appendix  C.  MATLAB  Post  Processing  Code 


C.l  CTH  Data  Extraction 

Output  data  from  CTH  are  stored  a  text  file  called  ’’hscth.”  The  output  file  is 
comma-delimited,  and  is  easily  opened  and  converted  to  tab-delimited  using  Excel. 
The  ’’hscth”  hie  includes  data  pertaining  to  the  CTH  cycle  number  and  current  step 
time.  This  is  found  in  the  second  and  third  column.  This  data  needs  to  be  removed 
before  passing  the  hie  through  the  MATLAB  post-processing  hie.  The  hrst  three  rows 
of  the  ’’hscth”  hie  needs  to  be  removed  as  well.  These  header  rows  give  the  titles  for 
each  column  and  are  unnecessary. 

After  the  two  columns  and  three  rows  are  removed,  the  data  set  should  consist 
of  columns  containg  data  in  this  order:  time,  x-position,  y-position,  von  Mises  stress, 
z-position,  xy-stress  deviator,  yy-stress  deviator,  xx-stress  deviator,  material  pressure 
(hydrostatic  stress),  volume  fraction  of  the  slipper,  plastic  strain  rate,  and  plastic 
strain  of  the  slipper.  The  default  filename  for  this  code  is  ’’cthData.txt”  but  this  can 
be  modified. 


C.2  MATLAB  Post  Processing  Code 

7.7,  CTH  DATA  POST  PROCESS  -  PLANE  STRAIN  EVALUATION 

X  Stephen  Meador  -  AFIT/GAE/ENY /10M-16 
4  7,  Master  ’  s  Student 

%  CTH  Slipper  Wear  Post  Process  Code 
%  Written  Sept  2009 -Mar  2010 

clear  all ;  close  all ;  clc 
9 

7.  HOW  TO  USE  THIS  POST-PROCESS  FILE: 

%{ 

This  file  is  divided  into  several  cells.  The  first  cell  clears  ... 
the 

workspace  and  closes  any  open  windows .  This  cell  also  has  two  . .  . 
variables 

14  that  must  be  defined  by  the  user:  "aspRad"  and  "velocity."  aspRad... 
is  the 

plane  strain  asperity  radius  with  units  of  micrometers , 
velocity  is  the 

collision  velocity  in  meters  per  second. 
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and 


The  second  cell  defines  the  poisson 1 s  ratio  of  the  slipper  ... 
material  and  the 

19  mesh  size  used  for  CTH  simulations  defined  with  units  of  ... 
centimeters 

squared.  This  code  assumes  a  uniform  CTH  cell  arrangement  where  . 
the  cells 

are  all  squares  of  equal  size.  This  is  important  for  calculating 
damage 

area  later  in  the  code.  Additionally,  the  directory  containing  .. 
the  CTH 

data  is  defined  based  on  the  simulation  asperity  radius  and  ... 
si ipper 

24  velocity  defined  in  the  first  cell. 

The  third  cell  imports  the  CTH  data,  and  the  fourth  cell  ... 
categorizes  as 

individual  arrays  and  matrices.  The  data  should  be  organized  such 
that 

each  row  represents  the  data  extracted  for  a  given  time  step  of  .. 
the  CTH 

29  simulation,  and  the  columns  are  the  data  extracted  from  the  CTH  .. 
tracer 

points.  The  user  should  note  the  order  in  which  the  variables  are 
arranged 

by  the  xxxLoc  variables  in  cell  four. 

Cell  five  calculates  the  sliding  distance  for  a  given  simulation. 
This 

34  distance  is  assumed  to  be  110%  of  the  asperity  radius  and  is  ... 
calculated 

based  on  the  velocity  and  simulation  time.  The  sixth  cell  ... 
calculates  the 

ZZ-deviatoric  stress  based  on  the  other  deviatoric  stresses  output 
from 

CTH.  The  seventh  cell  then  converts  all  stress  components  to  ... 
Pascals . 

39  The  eighth  cell  evaluates  the  strain  rates  at  every  tracer  point  . 
during 

the  simulation.  The  Johnson-Cook  constitutive  model  defines  the  . 
minimum 

strain  rate  as  0.002  1/s,  so  any  strain  rate  below  this  value  is  . 

reset  . 

Also,  if  any  strain  rate  exceeds  10~17  1/s  then  the  code  ... 
terminates 

because  the  stress  and  strain  curve  fits  have  not  been  evaluated  . 
for  data 

44  above  this  level. 

Cells  nine  and  ten  evaluate  the  stress  tensor  components  and  ... 
calculate  the 
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von  Mises  stresses,  respectively.  Cells  eleven  through,  fourteen 
e valuat  e 

the  various  failure  criteria.  And,  finally,  cell  fifteen  saves 
the  wear 

49  rate  data  to  text  files. 

aspRad  =  6;  "/,  microns 

54  velocity  =  700;  7  meters  per  second 

VMcr it  =  3 . 00  e9 ; 

Yo  =  1 . 897  e9  ; 

SigHEL  =  2 . 8664 e9  ; 

59 

tic 

77,  POISSON’S  RATIO,  MESH  SIZE, 

64  nu  =  0.283;  7  Poisson’s  ratio  of  material 

meshSize  =  1  .  Oe -4*  1 . 0  e -4 ;  7,  Area  of  a  single  mesh  cell  in 

~2 

if  velocity  <  100 

69  newDirectory  =  [’Data/00’  num2str (velocity)  ... 

’mps/0’  num2str ( aspRad )  ’micron’]; 
elseif  velocity  <1000 

newDirectory  =  [’Data/O’  num2str (velocity )  ... 

’mps/0’  num2str ( aspRad )  ’micron’]; 

74  else 

newDirectory  =  [’Data/’  num2str (velocity )  ... 

’mps/0’  num2str ( aspRad )  ’micron’]; 

end 

79  cd ( newDirectory ) 

di sp  (  ’  ’  ) 

7,7,  IMPORT  DATA 
84 

dataFile  =  ’cthData.txt’; 

data  =  load ( dat aFile ) ; 

89  disp(’Data  Imported  .  .  .  ’ ) 

77  CATEGORIZE  DATA 

time  =  data ( :  ,1)  ; 

94 

numCycles  =  length ( t ime ) ; 


cm . 
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numPoints  =  ( s ize ( dat a , 2) - 1 ) / 12 ; 


xPoints  =  zeros (numCycles .numPoints) ; 

99  yPoints  =  zeros (numCycles .numPoints) ; 

pressureData  =  zeros (numCycles , numPoints) 
vonMisesData  =  zeros (numCycles , numPoints) 
xxdevData  =  zeros (numCycles , numPoints ) ; 
yydevData  =  zeros (numCycles , numPoints ) ; 
104  xydevData  =  zeros (numCycles  , numPoints )  ; 
vfData  =  zeros (numCycles , numPoints ) ; 
srData  =  zeros (numCycles , numPoints ) ; 
strainData  =  zeros (numCycles , numPoints ) ; 
jcpData  =  zeros (numCycles , numPoints ) ; 

109 

xLoc  =  2 ; 
yLoc  =  3; 
vmLoc  =  5; 
xyLoc  =  6; 

114  yyLoc  =  7 ; 
xxLoc  =  8; 
pLoc  =  9 ; 
vfLoc  =  10; 
srLoc  =  11; 

119  sLoc  =  12; 
j  cpLoc  =  13  ; 

for  iter  =  1: numPoints 

xPoints (:, iter )  =  data ( : , xLoc ) ; 

124  yPoints (:, iter  )  =  data (:, yLoc )  ; 

pressureData (:, iter )  =  data ( : , pLoc ) ; 

vonMi sesDat a ( : , it er )  =  data (:, vmLoc ) ; 


xxdevDat  a ( : 

, iter ) 

=  data( 

, xxLoc )  ; 

yydevDat  a ( : 

,  iter ) 

=  data( 

, yyLoc)  ; 

129 

xydevDat  a ( : 

, iter ) 

=  data( 

, xyLoc )  ; 

vf Data (:, iter )  =  data (:, vfLoc ) ; 

srData (:, iter )  =  data (:, srLoc ) ; 

strainData (:, iter )  =  data ( : , sLoc ) ; 
j cpData ( : , iter )  =  data (:, j cpLoc ) ; 

134 


xLoc  = 

xLoc  + 

12; 

yLoc  = 

yLoc  + 

12; 

vmLoc 

=  vmLoc 

+  12; 

xyLoc 

=  xyLoc 

+  12; 

139 

yyLoc 

=  yyLoc 

+  12; 

xxLoc 

=  xxLoc 

+  12; 

pLoc  = 

pLoc  + 

12; 

vfLoc 

=  vfLoc 

+  12; 

srLoc 

=  srLoc 

+  12; 

144 

sLoc  = 

sLoc  + 

12; 

end 

j  cpLoc 

=  jcpLo 

c  +  12 
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disp (’Data  Categorized. . . ’  ) 

149 

ZZ  CALCULATE  DISTANCE  SLID 


distanceSlid  =  veloc ity *t ime ( end )* 1000 ;  Z  mm 

154  disp(’ Distance  Slid  Calculated...’) 

ZZ  CALCULATE  ZZDEV  (GIVEN  XXDEV ,  YYDEV ,  AND  POISSON > 
zzdevData  =  ( xxde vDat a+ yyde vDat a ) *nu ; 

159 

disp ( ’ ZZ  Deviator  Calculated...’) 

ZZ  CONVERT  DATA  TO  Pa 

164  pressureData  =  pressureData/10 ; 
xxdevData  =  xxdevData/10 ; 
yydevData  =  yydevData/10 ; 
xydevData  =  xydevData/10 ; 
zzdevData  =  zzdevData/10 ; 

169  jcpData  =  jcpData/10; 

vonMisesData  =  vonMisesData/10 ; 

ZZ  EVALUATE  STRAIN  RATES  FOR  ZEROS 


174  for  r  =  1  :  s  ize  (  srDat  a  ,  1 ) 

for  c  =  1 : size ( srData  ,  2) 
if  srData (r , c) <. 002 

srData(r,c)  =  .002; 


if  srData (r , c) >10el7 

disp (’Temp: ’) ,disp(temp) 
disp(’H  Vel :’), disp (velocity ) 
disp(’V  Vel :’), disp (vVel ) 

184  disp  (’ Row  :’),  disp  (r) 

dispC’Col: ’) , disp (c) 

disp(’Strain  Rate’) ,disp(srData(r,c)) 
error (’ Strain  Rate  Out  of  Range’) 

end 

189  end 

end 


ZZ  EVALUATE  STRAIN  AT  MAX  STRESS  FAILURE  AREA 

194 

failureSMS  =  zeros (numCycles , numPoints ) ; 
f ailureSumSMS  =  zeros (numCycles  ,  1)  ; 

199  A  =  2 . 24700e -2 ; 


RATIO) 
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B 

C 


-5 . 5160e -2 ; 
6 . 04400e -3 ; 


f ailureCr itSMS  =  A* ( srData . ~B)  +  C; 

204 

for  row=l : r 

for  col=l : c 

if  row>l  kk  f ai lur eSMS ( row - 1 , col ) =  =  1 
209  f ai lur eSMS ( row , col ) = 1 ; 


end 


214 


if  strainData (row , col) >=failureCritSMS (row , col) 
failureSMS (row , col) =1 ; 


end 

end 

219  end 

failureSMS  =  f ailureSMS . * vf Data ; 

for  iter  =  1 : length ( f ai lur eSumSMS ) 

224  f ai lur eSumSMS ( it er  ,  1 )  =  sum ( f ailureSMS ( iter  ,:))  ; 

end 

damAreaSMS  =  f ailureSumSMS*meshSize ; 

229  WR_SMS  =  100* damAreaSMS ( end )/ di st anc eSl id ; 

disp(’Strain  at  Max  Stress  Failure  Mechanism  Evaluated.  .  .  ’) 
Evaluate  CTH  J2P  data 

234 

f ailureSumVMS  =  zeros (numCycles  ,1)  ; 

VONMISESDATA  =  zer o s ( numCy c 1 e s , numPo int s ) ; 

239  for  row=l : r 

for  col=l : c 

if  row  > 1  kk  VONMISESDATA (row -1 , col ) ==1 
VONMISESDATA (row , col) =1 ; 

244 


end 

if  row>l  kk  vonMisesData (row -1 , col ) >=VMcrit 

249  VONMISESDATA (row , col) =1 ; 

end 


end 
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254 

259 

264 

269 

274 

279 

284 

289 

294 

299 


end 


VonMisesFracData  =  VONMISESDATA . * vf Data ; 

for  iter  =  1 : length ( f ai lureSumVMS ) 

f ai lur eSumVMS ( iter  , 1 )  =  sum ( VonMi sesFr acDat a ( iter  ,  :  ) )  ; 

end 

damAreaVMS  =  f ailureSumVMS*meshSize ; 

WR_VMS  =  100* damAreaVMS ( end )/ di st anc eSl id ; 

disp ( ’ VonMises  Stress  Failure  Mechanism  Evaluated.  .  .  ’) 

ZZ  SAVE  WEAR  RATES  TO  .  txt  FILE 
if  velocity  <  100 

fileName  =  [ 1 WearRates_00 ’  num2str (velocity )  ... 

’mps_0’  num2str ( aspRad )  ’micron.txt ’] ; 
elseif  velocity  <1000 

fileName  =  [ 1 WearRates_0 ’  num2str (velocity )  ... 

’mps_0’  num2str ( aspRad)  ’micron.txt ’] ; 

else 

fileName  =  [ ’ WearRates_ ’  num2str (velocity )  ... 

’mps_0’  num2str ( aspRad)  ’micron.txt’]; 

end 

fid=fopen(fileName , ’ wt ’) ; 

fprintf  (f  id  ,  ’  °/6 . 5e\t%6 .  5 e \t %6 . 5 e \t °/06 . 5 e \t  ’  ,  .  .  . 

WR_SMS *8 . 29 e -3  ,  .  .  . 

WR_ VMS *8 . 29e-3) ; 
fclose (fid)  ; 

disp(’Failure  Data  Saved...’) 

ZZ  END  PROGRAM 

disp (’ PROGRAM  COMPLETE...’) 

toe 

ZZ  Plot  Pressure  Time  Data  for  points  @  (x,y)  =  (0 .  0692 , 0 . 020f.) 

Pressure (1 : numCycles  ,1)  =  pressureData (1:101 ,424) -pressureData.  .  . 
(1 ,424)  ; 


f igur e 

plot (time , Pressure (1 : numCycles  ,  1 ) ) 
xlabel(’time  (s)’) 
ylabel (’ Pressure  (GPa)’) 

title (’ Pressure  at  (x,y)  =  (0.0692,0.0202)’) 
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grid  on 
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Appendix  D.  MATLAB  Code  for  Pressure  Along  a  Diagonal 

The  MATLAB  code  presented  in  this  appendix  is  used  to  plot  the  pressure  along  a 
30°  45°  and  60°  diagonal  from  the  point  of  contact  between  the  slipper  and  asperity 
with  respect  to  the  horizontal  rail.  The  CTH  input  deck  must  be  modified  before  the 
MATLAB  code  can  be  used.  The  tracer  input  set,  line  177  to  line  213  of  the  example 
CTH  input  deck  in  Appendix  B,  defines  the  initial  locations  of  the  data  points.  This 
section  must  be  modified  to  only  include  data  points  along  the  diagonals.  Replacing 
lines  177  to  213  with  the  following  lines  records  data  along  the  30°  diagonal.  This 
input  deck  must  be  run  three  times  to  record  the  data,  with  each  run  capturing  one 
diagonal.  The  asterisks  at  the  begining  of  a  line  comments  that  line  out  of  the  input 
deck.  To  capture  data  along  the  45°  diagonal,  an  asterisk  needs  to  be  added  to  line 
2  to  comment  the  30°  out,  and  the  asterisk  on  line  6  should  be  removed.  Removing 
the  asterisk  on  line  9  captures  the  data  along  the  60°  diagonal. 

D.l  Modified  Tracer  Input  Set 

tracer 

*  30  degrees 

add  693 . 4e -4 ,  200.6e-4  to  578.45e-4,  266.64e-4  n=500 
4 

*  45  degrees 

*  add  693.4e-4,  200.6e-4  to  600e-4,  294e-4  n=500 

*  60  degrees 

9  *  add  693 . 4e -4 ,  200.6e-4  to  627.36e-4,  315.53e-4  n=500 

endtracer 

When  the  three  simulations  are  finished,  the  first  two  columns  and  first  three 
rows  need  to  be  removed  in  Excel  using  the  process  outlined  in  Appendix  C.  The 
MATLAB  code  plots  the  change  in  pressure  along  the  three  diagonals  with  respect 
to  distance  for  each  time  step.  These  images  are  used  to  determine  the  clastic-plastic 
wave  speed  through  the  VascoMax  300  slipper  in  Section  4.4. 

D.2  MATLAB  Post  Processing  Code 

clear  all;  close  all;  clc 
7.7  VELOCITY  TO  PLOT 
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5 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 


veloc ity  =  500 ; 


angle  =  45; 

77.  LOAD  DATA 
if  velocity  <100 

datafile  =  [’00’  num2str (velocity)  ’ mps30degData . txt  ’  ]  ; 

elseif  velocity <1000 

datafile  =  [’O’  num2str (velocity )  ’ mps30degDat a . txt  ’  ]  ; 

else 

datafile  =  [num2str (velocity)  ’mps30degData.txt’]; 

end 

data30  =  load ( dataf i le ) ; 
if  velocity  <100 

datafile  =  [’00’  num2str (velocity)  ’mps45degData.txt’]; 

elseif  velocity <1000 

datafile  =  [’O’  num2str (velocity )  ’mps45degData.txt’]; 

else 

datafile  =  [num2str (velocity)  ’mps45degData.txt’]; 

end 

data45  =  load ( dataf i le ) ; 
if  velocity  <100 

datafile  =  [’00’  num2str (velocity)  ’mps60degData.txt’]; 

elseif  velocity <1000 

datafile  =  [’O’  num2str (velocity )  ’mps60degData.txt’]; 

else 

datafile  =  [num2str (velocity)  ’mps60degData.txt’]; 

end 

data60  =  load ( dataf i le ) ; 

7.7,  ORGANIZE  30  DEGREE  DATA 

time30  =  data30(:,l); 

nCycles30  =  size (data30 , 1) ; 
nPoints30  =  ( s ize ( dat a30 , 2 ) -1 ) /4 ; 

xPos30  =  2; 
yPos30  =  3; 
zPos30  =  4; 
pPos30  =  5; 

xData30  =  zeros (nCycles30 , nPoints30 ) ; 
yData30  =  zeros (nCycles30 , nPoints30 ) ; 
zData30  =  zeros (nCycles30 , nPoints30 ) ; 
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pData30  =  zeros (nCycles30 , nPoints30 ) ; 


for 

iter  =  1 

nPoints30 

xData30 ( 

, iter ) 

= 

dat a30  ( 

, xPos30 )  ; 

60 

yData30 ( 

, iter ) 

= 

dat a30  ( 

, yPos30 )  ; 

zData30 ( 

, iter ) 

= 

dat a30  ( 

, zPos30  )  ; 

pData30 ( 

, iter ) 

= 

dat a30  ( 

, pPos30  )  ; 

xPos30  = 

xPos30 

+ 

4; 

65 

yPos30  = 

yPos30 

+ 

4; 

zPos30  = 

zPos30 

+ 

4; 

pPos30  = 

pPos30 

+ 

4; 

end 


70  Z7,  ORGANIZE  45  DEGREE  DATA 

time45  =  dat a45 ( :  ,  1 )  ; 

nCycles45  =  size (data45 , 1) ; 

75  nPoints45  =  ( s ize ( data45  ,  2 )  -1 ) /4 ; 


xPos45  =  2 
yPos45  =  3 
zPos45  =  4 
80  pPos45  =  5 


xData45  =  zeros (nCycles45 , nPoints45 ) ; 
yData45  =  zeros (nCycles45 , nPoints45 ) ; 
zData45  =  zeros (nCycles45 , nPoints45 ) ; 
85  pData45  =  zeros (nCycles45 , nPoints45 ) ; 


f  or 

iter  =  1 

: nPoints45 

xData45 ( 

: , iter ) 

= 

dat a45  ( 

, xPos45 )  ; 

yData45 ( 

: , iter ) 

= 

dat a45  ( 

, yPos45 )  ; 

90 

zData45 ( 

: , iter ) 

= 

dat a45  ( 

, zPos45 )  ; 

pData45 ( 

: , iter ) 

= 

dat a45  ( 

, pPos45 )  ; 

xPos45  = 

xPos45 

+ 

4; 

yPos45  = 

yPos45 

+ 

4; 

95 

zPos45  = 

zPos45 

+ 

4; 

pPos45  = 

pPos45 

+ 

4; 

end 

ZZ 

ORGANIZE 

60  DEGREE 

DATA 

100 


time60  =  data60(:,l); 


nCycles60  =  size (data60 , 1) ; 
nPoints60  =  ( size ( data60  ,  2)  -1) /4 ; 

105 

xPos60  =  2; 
yPos60  =  3; 
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zPos60  =  4; 
pPos60  =  5; 

110 


xData60 

yData60 

zData60 

pData60 


zeros (nCycles60 
zeros (nCycles60 
zeros (nCycles60 
zeros (nCycles60 


nPoints60 )  ; 
nPoints60 )  ; 
nPoints60  )  ; 
nPoints60  )  ; 
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for  iter  =  1 

nPoints60 

xData60 ( 

, iter ) 

= 

data60  ( 

, xPos60 )  ; 

yData60 ( 

, iter ) 

= 

data60  ( 

, yPos60 )  ; 

zData60 ( 

, iter ) 

= 

data60  ( 

, zPos60  )  ; 

120  pData60 ( 

, iter ) 

= 

data60  ( 

, pPos60  )  ; 

xPos60  = 

xPos60 

+ 

4; 

yPos60  = 

yPos60 

+ 

4; 

zPos60  = 

zPos60 

+ 

4; 

125  pPos60  = 

pPos60 

+ 

4; 

end 


y°/,  CONVERT  UNITS 


130  plnt30 

pDat  a30 

pDat  a30 
135  pData30 

xDat  a30 
yDat  a30 


mean(pData30 (1  ,  : ) )  ; 

pData30 -plnt30 ; 

pDat  a30 / 1 0 ; 
pDat  a30 / 10 “ 9 ; 

xDat  a30  *  10 " 4 ; 
yDat  a30  *  10 " 4 ; 


140  plnt45 

pDat  a45 

pDat  a45 
145  pData45 

xDat  a45 
yDat  a45 


mean(pData45 (1  ,  : ) )  ; 

pData45 -plnt45 ; 

pDat  a45 / 1 0 ; 
pData45 /10~9 ; 

xDat  a45  *  10  ~  4 ; 
yData45*10~4; 


150  plnt60 

pData60 

pData60 
155  pData60 

xData60 

yData60 


mean(pData60 (1  ,  : ) )  ; 

pData60 -plnt60 ; 

pDat  a60 / 1 0 ; 
pData60 / 1 0  ~  9 ; 

xDat  a60  *  10 " 4 ; 
yDat  a60  *  10 " 4 ; 


%  Initial  Pressure 

/  Pressure  Change 

%  Pa 
%  GPa 

°/o  microns 
"/,  microns 

°/o  Initial  Pressure 

%  Pressure  Change 

/  Pa 
y  GPa 

°/,  microns 
y  microns 

y  Initial  Pressure 

°/o  Pressure  Change 

°/o  Pa 
/o  GPa 

y  microns 
y  microns 
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160 

165 

170 

175 

180 

185 

190 

195 

200 

205 

210 


XX  DEFINE  LENGTH 


minX30  = 
maxX30  = 
minY30  = 
maxY30  = 

deltaX30 

deltaY30 

length30 

tracerPo 

minX45  = 
maxX45  = 
minY45  = 
maxY45  = 

deltaX45 

deltaY45 

length45 

tracerPo 

minX60  = 
maxX60  = 
minY60  = 
maxY60  = 

deltaX60 

deltaY60 

length60 

tracerPo 

XX  PLOT 


min ( xDat  a30 ( 1  , 

:))  ; 

max ( xDat  a30 ( 1 , 

:))  ; 

min ( yDat  a30 ( 1 , 

:))  ; 

max ( yData30 ( 1 , 

:)  )  ; 

=  maxX30  -  minX30 ; 

=  maxY30  -  minY30  ; 

=  sqrt ( deltaX30 ~2  +  deltaY30~2); 

ints30  =  linspace (0 , length30 , nPoints30 ) ; 


min ( xDat  a45 (1  , 

:))  ; 

max ( xDat  a45 (1 , 

:))  ; 

min ( yData45 ( 1 , 

:))  ; 

max (yData45 (1 , 

:)  )  ; 

=  maxX45  -  minX45 ; 

=  maxY45  -  minY45  ; 

=  sqrt ( deltaX45 ~2  +  deltaY45~2); 

ints45  =  linspace  (0  ,  length.45  ,  nPoints45  )  ; 


min ( xData60 ( 1 , 

:)  )  ; 

max ( xData60 ( 1 , 

:)  )  ; 

min ( yData60 ( 1 , 

:)  )  ; 

max ( yData60 ( 1 , 

:)  )  ; 

=  maxX60  -  minX60 ; 

=  maxY60  -  minY60; 

=  sqrt ( deltaX60 ~2  +  deltaY60~2); 

ints60  =  linspace (0 , length60 , nPoints60 ) ; 


figNum  =  0; 

plotTracers  =  0; 
if  plotTracers 

for  index  =  1 X : 1 00 : nCy cl es 

x  =  [500  693.4] ; 

y30  =  [200*tand (30) +200 . 6  200.6] 
y45  =  [200*tand (45) +200 . 6  200.6] 
y60  =  [200*tand (60) +200 . 6  200.6] 
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testX 

=  xDat a30 ( index  , 

:)  ; 

testY 

=  yDat a30 ( index  , 

:)  ; 

215 

testP 

=  pDat a30 ( index  , 

:)  ; 

figNum  =  figNum  +  1; 
f igur e ( f igNum  ) 
hold  on 

220  plot (testX , testY  ,  ’  d  ’  ) 

plot (x , y30  ,  ’  k  ’  ) 
plot (x , y45  ,  ’  k  ’  ) 
plot (x , y60  ,  ’  k  ’  ) 
axis ( ’ square  ’  ) 

225  xLim (  [500  700]) 

yLim (  [200  400]) 

end 

end 

230 

plotPressures  =  1; 
if  plotPressures 

for  index  =  1 : 1 : min ( [nCycles30  nCycles45  nCycles60]) 

235  titleText  =  sprintf  (  1  °/„4 .  Of  m/s  I  Pressure  Change  Along  ... 

Diagonal  at  '/,  1 . 2  e  seconds  1  ,  velocity  ,  time30  (  index)  )  ; 

p30  =  pDat a30 ( index  ,:)  ; 
p45  =  pDat a45 ( index  ,:)  ; 
p60  =  pData60 ( index  ,:)  ; 

240 

figNum  =  figNum  +  1; 
f igur e ( f igNum  ) 
subplot ( 1 , 7  ,  [3  7]) 
hold  on 

245  plot(tracerPoints30  ,p30,  ’  r  ’  ,  ’LineWidth’  ,2) 

plot(tracerPoints45  ,p45,  ’  k  ’  ,  ’LineWidth1  ,2) 
plot(tracerPoints60  ,p60,  ’  b  ’  ,  ’LineWidth1  ,2) 

%  plot  ( [6  6],  [-10  15],  >i - LineWidth  ’  ,2) 

xlim ( [min (min ( [tracerPoints30 ;  tracerPoints45 ;  ... 

tracerPoints60 ] ) )  max (max ( [tracerPoints30 ;  ... 

tracerPoints45 ;  tracerPoints60] ) ) ] ) 

250  ylim  (  [-15  15]) 

"/,  title  ({ [num2str  (velocity )  ’  m/s  I  Pressure  Change  Along 

Diagonal  at  Time  ’  num2str (time30 (index) )  ’  seconds  ’]}) 

title (titleText) 

xlabel (’ Distance  Along  Diagonal  (\mum) ’) 
ylabel ( ’ \DeltaP  (GPa)’) 

255  grid  on 

legend(’30  deg’, ’45  deg ’ , ’ 60  deg ’,’ Location ’,’ NorthEast ’ ) 
subplot (1  ,7,1) 

bar (( index -1) / (min ( [nCycles30  nCycles45  nCycles60] ) -1) . .  . 
*100 ,  ’k ’  ,  ’ BarWidth ’  ,1) 
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ylim (  [0  100] ) 

260  ylabel ( ’ Percentage  of  Simulation  Time  Completed’) 

if  index<10 

if  velocity  <100 

saveas (gcf  ,  [  ’  00  ’  num2str (velocity )  ’mpsOO’  num2str 
( index )  ’  . bmp  ’  ] ) 

265  elseif  velocity <1000 

saveas (gcf  ,[’ 0  ’  num2str (velocity)  ’mpsOO’  num2str( 
index )  ’  . bmp  ’  ]  ) 

else 

saveas (gcf  ,[’ 0  ’  num2str (velocity )  ’mpsOO’  num2str( 
index )  ’  . bmp  ’  ] ) 

end 

270  elseif  index<100 

if  velocity  <100 

saveas (gcf  ,[’  00  ’  num2str (velocity )  ’mpsO’  num2str( 
index )  ’  . bmp  ’  ]  ) 

elseif  velocity <1000 

saveas (gcf  ,[’  0  ’  num2str (velocity )  ’mpsO’  num2str(. 
index )  ’  . bmp  ’  ] ) 

275  else 

saveas (gcf  ,[’ 0  ’  num2str (velocity )  ’mpsO’  num2str(. 
index )  ’  . bmp  ’  ] ) 

end 

else 

if  velocity  <100 

280  saveas (gcf ,[’ 00 ’  num2str (velocity )  ’mps ’  num2str(. 

index )  ’  . bmp  ’  ]  ) 

elseif  velocity <1000 

saveas (gcf  ,[’  0  ’  num2str (velocity )  ’mps’  num2str(.. 
index )  ’  . bmp  ’  ] ) 

else 

saveas (gcf  ,[’ 0  ’  num2str (velocity )  ’mps’  num2str(.. 
index )  ’  . bmp  ’  ] ) 

285  end 

end 

close ( gcf ) 


290  end 

end 

figNum  =  figNum  +  1; 
figure ( f igNum ) 

295  plot(time30 ,pData30(:  ,114)  ,  ’k’  ,  ’LineWidth’  ,2) 
xlim([0  max ( t ime30 ) ] ) 
ylim (  [-15  15]) 

title (’ Pressure  at  30  \mum  Along  30~{\circ}  Diagonal’) 
xlabel(’Time  (s)’) 

300  ylabel (’ \DeltaP  (GPa)’) 


112 


figNum  =  figNum  +  1; 
figure (figNum) 

plot(time30  ,pData30(:  ,227)  ,  '  k  ’  ,  'LineWidth'  ,2) 

305  xlim([0  max ( t ime30 ) ] ) 
ylim (  [-15  15]) 

title (’ Pressure  at  60  \mum  Along  30~{\circ}  Diagonal') 
xlabel(’Time  (s)') 
ylabel (' \DeltaP  (GPa)’) 

310 

figNum  =  figNum  +  1; 
figure (figNum) 

plot(time45  ,pData45(:  ,115)  ,  '  k  ’  ,  'LineWidth'  ,2) 
xlim([0  max ( t ime45 ) ] ) 

315  ylim ( [-15  15]  ) 

title (’ Pressure  at  30  \mum  Along  45~{\circ}  Diagonal') 
xlabel('Time  (s)') 
ylabel ( ' \DeltaP  (GPa)’) 

320  figNum  =  figNum  +  1; 
figure (figNum) 

plot(time45  ,pData45(:  ,229)  ,  'k'  ,  'LineWidth'  ,2) 
xlim([0  max ( t ime45 ) ] ) 
ylim (  [-15  15]) 

325  title (’ Pressure  at  60  \mum  Along  45~{\circ}  Diagonal') 
xlabel('Time  (s)') 
ylabel (' \DeltaP  (GPa)’) 

figNum  =  figNum  +  1; 

330  f  igur  e  (  f  igNum  ) 

plot(time60  ,pData60(:  ,114)  ,  'k'  ,  'LineWidth.'  ,2) 
xlim([0  max ( t ime60 ) ] ) 
ylim (  [-15  15]) 

title (’ Pressure  at  30  \mum  Along  60~{\circ}  Diagonal') 
335  xlabel('Time  (s)') 

ylabel (' \DeltaP  (GPa)’) 

figNum  =  figNum  +  1; 
figure (figNum) 

340  plot(time60 ,pData60(:  ,227)  ,  'k'  ,  'LineWidth'  ,2) 
xlim([0  max ( t ime60 ) ] ) 
ylim (  [-15  15]) 

title (’ Pressure  at  60  \mum  Along  60~{\circ}  Diagonal') 
xlabel('Time  (s)') 

345  ylabel (' \DeltaP  (GPa)’) 
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